library(raster)

library(plyr)
library(dplyr)

library(ggplot2)
library(viridisLite)
library(colorspace)
library(RColorBrewer)

library(sf)
library(stringr)

library(ggpubr)
library(kableExtra)

memory.limit(size = 160000)
## [1] Inf

#display.brewer.pal(n = 11, name ="RdYlGn")
# colorvec <- brewer.pal(n = 11, name ="RdYlGn")
# col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

#handmade
col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")

1 Data

Land-sea mask and IPCC region are taken from Iturbide et al. 2020: https://github.com/SantanderMetGroup/ATLAS/tree/devel/reference-grids

Elevation data is taken from Worldclim: https://www.worldclim.org/data/worldclim21.html

land_mask <- raster("Masks/land_sea_mask_1degree.nc4") 
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land


elev <- raster("Worldclim/wc2.1_10m/wc2.1_10m_elev.tif") 
elev.df <- elev %>% projectRaster(to = land_mask) %>% as.data.frame(xy = T) %>% setNames(c("lon","lat", "z"))

ipcc_regions <- shapefile("Masks/IPCC-WGI-reference-regions-v4.shp") %>% spTransform(crs("EPSG:4326")) 
ipcc_regions.raster <- ipcc_regions %>% rasterize(land_mask)
ipcc_regions.df <- as.data.frame(ipcc_regions.raster, xy = T) %>% setNames(c("lon","lat","Continent","Type","Name","Acronym"))

CMIP6 data is dowloaded from https://aims2.llnl.gov/search

The data was downloaded by model and extracted by model, period and SSP with the “data-prep-xxx.R” scripts, then grouped by variable (“bind-tables_xxx.R”), on GRICAD.

tas.all_annual <- read.table("CMIP6/tas/tas.all_annual.txt")
pr.all_annual <- read.table("CMIP6/pr/pr.all_annual.txt")
sfcWind.all_annual <- read.table("CMIP6/sfcWind/sfcWind.all_annual.txt")
hfls.all_annual <- read.table("CMIP6/hfls/hfls.all_annual.txt")
hfss.all_annual <- read.table("CMIP6/hfss/hfss.all_annual.txt")
hurs.all_annual <- read.table("CMIP6/hurs/hurs.all_annual.txt")

cmip6_annual <- tas.all_annual %>% merge(pr.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_annual, "cmip6_annual.txt")
tas.all_january <- read.table("CMIP6/tas/tas.all_january.txt")
pr.all_january <- read.table("CMIP6/pr/pr.all_january.txt")
sfcWind.all_january <- read.table("CMIP6/sfcWind/sfcWind.all_january.txt")
hfls.all_january <- read.table("CMIP6/hfls/hfls.all_january.txt")
hfss.all_january <- read.table("CMIP6/hfss/hfss.all_january.txt")
hurs.all_january <- read.table("CMIP6/hurs/hurs.all_january.txt")

cmip6_january <- tas.all_january %>% merge(pr.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_january, "cmip6_january.txt")
tas.all_february <- read.table("CMIP6/tas/tas.all_february.txt")
pr.all_february <- read.table("CMIP6/pr/pr.all_february.txt")
sfcWind.all_february <- read.table("CMIP6/sfcWind/sfcWind.all_february.txt")
hfls.all_february <- read.table("CMIP6/hfls/hfls.all_february.txt")
hfss.all_february <- read.table("CMIP6/hfss/hfss.all_february.txt")
hurs.all_february <- read.table("CMIP6/hurs/hurs.all_february.txt")

cmip6_february <- tas.all_february %>% merge(pr.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_february, "cmip6_february.txt")
tas.all_march <- read.table("CMIP6/tas/tas.all_march.txt")
pr.all_march <- read.table("CMIP6/pr/pr.all_march.txt")
sfcWind.all_march <- read.table("CMIP6/sfcWind/sfcWind.all_march.txt")
hfls.all_march <- read.table("CMIP6/hfls/hfls.all_march.txt")
hfss.all_march <- read.table("CMIP6/hfss/hfss.all_march.txt")
hurs.all_march <- read.table("CMIP6/hurs/hurs.all_march.txt")

cmip6_march <- tas.all_march %>% merge(pr.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_march, "cmip6_march.txt")
tas.all_april <- read.table("CMIP6/tas/tas.all_april.txt")
pr.all_april <- read.table("CMIP6/pr/pr.all_april.txt")
sfcWind.all_april <- read.table("CMIP6/sfcWind/sfcWind.all_april.txt")
hfls.all_april <- read.table("CMIP6/hfls/hfls.all_april.txt")
hfss.all_april <- read.table("CMIP6/hfss/hfss.all_april.txt")
hurs.all_april <- read.table("CMIP6/hurs/hurs.all_april.txt")

cmip6_april <- tas.all_april %>% merge(pr.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_april, "cmip6_april.txt")
tas.all_may <- read.table("CMIP6/tas/tas.all_may.txt")
pr.all_may <- read.table("CMIP6/pr/pr.all_may.txt")
sfcWind.all_may <- read.table("CMIP6/sfcWind/sfcWind.all_may.txt")
hfls.all_may <- read.table("CMIP6/hfls/hfls.all_may.txt")
hfss.all_may <- read.table("CMIP6/hfss/hfss.all_may.txt")
hurs.all_may <- read.table("CMIP6/hurs/hurs.all_may.txt")

cmip6_may <- tas.all_may %>% merge(pr.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_may, "cmip6_may.txt")
tas.all_june <- read.table("CMIP6/tas/tas.all_june.txt")
pr.all_june <- read.table("CMIP6/pr/pr.all_june.txt")
sfcWind.all_june <- read.table("CMIP6/sfcWind/sfcWind.all_june.txt")
hfls.all_june <- read.table("CMIP6/hfls/hfls.all_june.txt")
hfss.all_june <- read.table("CMIP6/hfss/hfss.all_june.txt")
hurs.all_june <- read.table("CMIP6/hurs/hurs.all_june.txt")

cmip6_june <- tas.all_june %>% merge(pr.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_june, "cmip6_june.txt")
tas.all_july <- read.table("CMIP6/tas/tas.all_july.txt")
pr.all_july <- read.table("CMIP6/pr/pr.all_july.txt")
sfcWind.all_july <- read.table("CMIP6/sfcWind/sfcWind.all_july.txt")
hfls.all_july <- read.table("CMIP6/hfls/hfls.all_july.txt")
hfss.all_july <- read.table("CMIP6/hfss/hfss.all_july.txt")
hurs.all_july <- read.table("CMIP6/hurs/hurs.all_july.txt")

cmip6_july <- tas.all_july %>% merge(pr.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_july, "cmip6_july.txt")
tas.all_august <- read.table("CMIP6/tas/tas.all_august.txt")
pr.all_august <- read.table("CMIP6/pr/pr.all_august.txt")
sfcWind.all_august <- read.table("CMIP6/sfcWind/sfcWind.all_august.txt")
hfls.all_august <- read.table("CMIP6/hfls/hfls.all_august.txt")
hfss.all_august <- read.table("CMIP6/hfss/hfss.all_august.txt")
hurs.all_august <- read.table("CMIP6/hurs/hurs.all_august.txt")

cmip6_august <- tas.all_august %>% merge(pr.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_august, "cmip6_august.txt")
tas.all_september <- read.table("CMIP6/tas/tas.all_september.txt")
pr.all_september <- read.table("CMIP6/pr/pr.all_september.txt")
sfcWind.all_september <- read.table("CMIP6/sfcWind/sfcWind.all_september.txt")
hfls.all_september <- read.table("CMIP6/hfls/hfls.all_september.txt")
hfss.all_september <- read.table("CMIP6/hfss/hfss.all_september.txt")
hurs.all_september <- read.table("CMIP6/hurs/hurs.all_september.txt")

cmip6_september <- tas.all_september %>% merge(pr.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_september, "cmip6_september.txt")
tas.all_october <- read.table("CMIP6/tas/tas.all_october.txt")
pr.all_october <- read.table("CMIP6/pr/pr.all_october.txt")
sfcWind.all_october <- read.table("CMIP6/sfcWind/sfcWind.all_october.txt")
hfls.all_october <- read.table("CMIP6/hfls/hfls.all_october.txt")
hfss.all_october <- read.table("CMIP6/hfss/hfss.all_october.txt")
hurs.all_october <- read.table("CMIP6/hurs/hurs.all_october.txt")

cmip6_october <- tas.all_october %>% merge(pr.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_october, "cmip6_october.txt")
tas.all_november <- read.table("CMIP6/tas/tas.all_november.txt")
pr.all_november <- read.table("CMIP6/pr/pr.all_november.txt")
sfcWind.all_november <- read.table("CMIP6/sfcWind/sfcWind.all_november.txt")
hfls.all_november <- read.table("CMIP6/hfls/hfls.all_november.txt")
hfss.all_november <- read.table("CMIP6/hfss/hfss.all_november.txt")
hurs.all_november <- read.table("CMIP6/hurs/hurs.all_november.txt")

cmip6_november <- tas.all_november %>% merge(pr.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_november, "cmip6_november.txt")
tas.all_december <- read.table("CMIP6/tas/tas.all_december.txt")
pr.all_december <- read.table("CMIP6/pr/pr.all_december.txt")
sfcWind.all_december <- read.table("CMIP6/sfcWind/sfcWind.all_december.txt")
hfls.all_december <- read.table("CMIP6/hfls/hfls.all_december.txt")
hfss.all_december <- read.table("CMIP6/hfss/hfss.all_december.txt")
hurs.all_december <- read.table("CMIP6/hurs/hurs.all_december.txt")

cmip6_december <- tas.all_december %>% merge(pr.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_december, "cmip6_december.txt")

2 Aridity Index

2.1 De Martonne

cmip6 <- read.table("cmip6.ref.txt") # created in "figures_article.Rmd", already converted in mm and °C

breaks.martonne <- c(0,5,10,20,"30","40",Inf)
cat.martonne <- c("[0,5]" = "Desert", "(5,10]"= "Arid", "(10,20]" = "Semi-arid","(20,30]" = "Temperate", "(30,40]" = "Humid", "(40,Inf]" = "Forest")
col.martonne <- c("Desert" =  "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Temperate" = "#ABD9E9", "Humid"= "#4575B4","Forest"=  "#DAE2F7")

cmip6$AIm <- with(cmip6,spr.ref/(t.ref)) # pr in mm/y, tas in Ceslsius

cmip6$cat.AIm <- cmip6$AIm %>% cut(breaks.martonne, include.lowest = T) %>% revalue(cat.martonne)

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1970_2000" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AIm))+
    scale_fill_manual(values = col.martonne, na.translate = F)+
    labs(title = paste(i, "1970-2000", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

2.2 AI FAO (with Penman-Monteith)

ea from Rh:

ea = 6.108(RH/100)exp(17.27*T/(T+237.3))

Refs: https://www.nature.com/articles/s41598-023-40499-6 : Allen and Tetens

All variables are calculated by month, so that we know the monthly evapotranspiration. The monthly values are summed to obtain the annual values.

id_vars <- c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")

cmip6_annual <- mutate(read.table("cmip6_annual.txt"),
                       t = tas - 273.15,
                       Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                       P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                       es = 0.61078*exp(17.27*t/(t+237.3)),
                       ea = (hurs/100)*es,
                       delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                       gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                       ET0_annual = ((0.408*delta*Rn+gamma*900/((t)+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind)))) 

cmip6_01 <- mutate(read.table("cmip6_january.txt"),
                   t = tas - 273.15,
                   pr_01 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_01 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind)))) # Evapotranspiration, G considered to be negligible. in mm/day. Wind speed at 10 m is converted to wind speed at 2 m by multiplying by 0.748

cmip6_02 <- mutate(read.table("cmip6_february.txt"),
                   t = tas - 273.15,
                   pr_02 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_02 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_03 <- mutate(read.table("cmip6_march.txt"),
                   t = tas - 273.15,
                   pr_03 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_03 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_04 <- mutate(read.table("cmip6_april.txt"),
                   t = tas - 273.15,
                   pr_04 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_04 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_05 <- mutate(read.table("cmip6_may.txt"),
                   t = tas - 273.15,
                   pr_05 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_05 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_06 <- mutate(read.table("cmip6_june.txt"),
                   t = tas - 273.15,
                   pr_06 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_06 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_07 <- mutate(read.table("cmip6_july.txt"),
                   t = tas - 273.15,
                   pr_07 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_07 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_08 <- mutate(read.table("cmip6_august.txt"),
                   t = tas - 273.15,
                   pr_08 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_08 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_09 <- mutate(read.table("cmip6_september.txt"),
                   t = tas - 273.15,
                   pr_09 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_09 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_10 <- mutate(read.table("cmip6_october.txt"),
                   t = tas - 273.15,
                   pr_10 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_10 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_11 <- mutate(read.table("cmip6_november.txt"),
                   t = tas - 273.15,
                   pr_11 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_11 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_12 <- mutate(read.table("cmip6_december.txt"),
                   t = tas - 273.15,
                   pr_12 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_12 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

et0_monthly <- select(cmip6_annual, id_vars) %>%
  merge(select(cmip6_01, c(id_vars, "ET0_01", "pr_01")), by = id_vars) %>%
  merge(select(cmip6_02, c(id_vars, "ET0_02", "pr_02")), by = id_vars) %>%
  merge(select(cmip6_03, c(id_vars, "ET0_03", "pr_03")), by = id_vars) %>%
  merge(select(cmip6_04, c(id_vars, "ET0_04", "pr_04")), by = id_vars) %>%
  merge(select(cmip6_05, c(id_vars, "ET0_05", "pr_05")), by = id_vars) %>%
  merge(select(cmip6_06, c(id_vars, "ET0_06", "pr_06")), by = id_vars) %>%
  merge(select(cmip6_07, c(id_vars, "ET0_07", "pr_07")), by = id_vars) %>%
  merge(select(cmip6_08, c(id_vars, "ET0_08", "pr_08")), by = id_vars) %>%
  merge(select(cmip6_09, c(id_vars, "ET0_09", "pr_09")), by = id_vars) %>%
  merge(select(cmip6_10, c(id_vars, "ET0_10", "pr_10")), by = id_vars) %>%   
  merge(select(cmip6_11, c(id_vars, "ET0_11", "pr_11")), by = id_vars) %>%
  merge(select(cmip6_12, c(id_vars, "ET0_12", "pr_12")), by = id_vars)

et0_monthly$spr <- with(et0_monthly, pr_01*31 + pr_02 * 28 + pr_03 * 31 + pr_04 * 30 + pr_05 * 31 + pr_06 * 30 +
                          pr_07 * 31 + pr_08 * 31 + pr_09 * 30 + pr_10 * 31 + pr_11 * 30 + pr_12 * 31)
et0_monthly$sET0 <- with(et0_monthly, ET0_01*31 + ET0_02 * 28 + ET0_03 * 31 + ET0_04 * 30 + ET0_05 * 31 + ET0_06 * 30 +
                           ET0_07 * 31 + ET0_08 * 31 + ET0_09 * 30 + ET0_10 * 31 + ET0_11 * 30 + ET0_12 * 31)

cmip6 <- merge(cmip6_annual, et0_monthly, by = id_vars)
cmip6$AI <- with(cmip6, abs(spr/sET0))  
cmip6$AI_annual <- with(cmip6, abs((pr*60*60*24*365)/(ET0_annual*365)))

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
#col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

cmip6$cat.AI <- cmip6$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6$cat.AI[which(cmip6$sET0 < 400)] <- "Cold"
cmip6$cat.AI %>% unique()

cmip6$cat.AI <- cmip6$AI %>% cut(breaks = breaks.unesco) %>% 
  revalue(cat.unesco)


cmip6$cat.AI_annual <- cmip6$AI_annual %>% cut(breaks = breaks.unesco) %>% 
  revalue(cat.unesco)

cmip6$cat.AI_annual[which(cmip6$ET0_annual*365 < 400)] <- "Cold"

write.table(cmip6, "cmip6_AI.txt")

3 Maps AI by model, all % calculated without Antarctica

cmip6 <- read.table("cmip6_AI.txt")

breaks.unesco <- c(-Inf,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(-Inf,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid")
cmip6$cat.AI <- factor(cmip6$cat.AI, levels =  c("Hyperarid","Arid","Semi-arid","Dry subhumid", "Humid","Cold"))

3.1 1850-1880

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1850_1880" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "1850-1880", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")


sdmm <- cmip6 %>% subset(period == "1850_1880"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))


tab <- cmip6 %>% subset(period == "1850_1880"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)


knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 5.1 9.0 13.2 6.1 33.4 38.0 28.3 0.3
CESM 5.2 10.8 14.7 5.0 35.7 35.6 28.5 0.3
CMCC 5.6 7.2 9.5 5.3 27.6 42.1 30.1 0.3
CMCC-ESM2 6.0 8.2 7.5 2.8 24.5 24.0 51.2 0.3
CNRM 5.7 12.1 11.4 6.6 35.8 38.0 25.9 0.2
EC-Earth3 8.6 9.0 11.0 4.5 33.1 36.7 30.1 0.2
FGOALS 5.6 12.8 15.9 7.2 41.5 29.9 28.4 0.2
GFDL-ESM4 5.6 9.3 8.7 4.1 27.7 37.0 34.9 0.2
INM 4.4 8.2 13.4 5.8 31.8 40.6 27.4 0.3
INM-CM5 3.5 8.8 11.5 5.6 29.4 43.2 27.1 0.3
MPI 8.5 10.9 10.1 4.4 33.9 34.5 31.4 0.2
MRI 6.3 10.6 8.7 3.8 29.4 39.6 30.9 0.3
NorESM-2-MM 5.1 10.4 18.4 5.9 39.8 35.4 24.4 0.3
mean 5.8 9.8 11.8 5.2 32.6 36.5 30.7 0.3
sd 1.4 1.6 3.2 1.2 7.4 5.1 6.7 0.1

3.2 1970-2000

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1970_2000" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "1970-2000", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

library(ggrepel) 

pie <- cmip6 %>% subset(period == "1970_2000" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  select(c("source", "cat.AI")) %>%
  reshape2::melt(id.vars = c("source")) %>%
  group_by(source, value) %>%
  summarise(count = n()) %>%
  mutate(perc = round(count/sum(count)*100, 1),
         lab.y = rev(cumsum(rev(count))) - count*0.5) %>% 
  ungroup() 

pie_list <- list()

for(i in unique(pie$source)){
  p <- ggplot(subset(pie, source == i))+
    geom_bar(aes(x = "", y = count, fill = value), width = 1, stat = "identity")+
    coord_polar("y", start = 0)+
    scale_fill_manual(values = col.cat)+
    geom_label_repel(aes(x = "", y = lab.y, label = perc), nudge_x = 0.5)+
    labs(fill = "", title = i)+
    theme_void()+theme(legend.position = "right")
  pie_list[[i]] <- p
}

ggpubr::ggarrange(plotlist = pie_list, nrow = 4, ncol = 4, common.legend = T, legend = "bottom")

pie <- cmip6 %>% subset(period == "1970_2000" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  select(c("cat.AI")) %>%
  table() %>% as.data.frame() %>%
  mutate(perc = round(Freq/sum(Freq)*100, 1),
   lab.y = rev(cumsum(rev(Freq))) - Freq*0.5) 

p <- ggplot(pie)+
    geom_bar(aes(x = "", y = Freq, fill = cat.AI), width = 1, stat = "identity")+
    coord_polar("y", start = 0)+
    scale_fill_manual(values = col.cat)+
    geom_label_repel(aes(x = "", y = lab.y, label = perc), nudge_x = 0.5)+
    labs(fill = "", title = "Multimodel CMIP6")+
    theme_void()+theme(legend.position = "right")

print(p)


sdmm <- cmip6 %>% subset(period == "1970_2000"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "1970_2000"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)


knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 5.2 8.6 13.2 6.3 33.3 37.4 29.1 0.3
CESM 5.3 10.7 14.1 5.2 35.3 35.0 29.5 0.3
CMCC 5.1 7.4 10.1 5.4 28.0 42.2 29.5 0.3
CMCC-ESM2 6.4 7.9 7.7 2.5 24.5 24.4 50.8 0.3
CNRM 5.8 12.2 11.3 6.5 35.8 38.3 25.7 0.2
EC-Earth3 8.4 9.5 10.8 4.4 33.1 35.7 31.0 0.2
FGOALS 5.9 10.9 10.0 4.9 31.7 25.1 42.9 0.2
GFDL-ESM4 5.8 8.9 8.9 4.1 27.7 35.3 36.8 0.2
INM 4.4 8.4 13.3 5.7 31.8 40.1 27.8 0.3
INM-CM5 3.4 8.5 12.0 5.6 29.5 43.0 27.3 0.3
MPI 8.7 11.1 10.2 4.3 34.3 34.0 31.5 0.2
MRI 6.1 10.8 9.2 3.9 30.0 38.9 30.9 0.3
NorESM-2-MM 5.2 10.1 17.1 5.6 38.0 33.4 28.3 0.3
mean 5.8 9.6 11.4 5.0 31.8 35.6 32.4 0.3
sd 1.4 1.5 2.5 1.1 6.5 5.7 7.1 0.1

3.3 1985-2010

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1985_2015" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "1985-2010", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")


sdmm <- cmip6 %>% subset(period == "1985_2015"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))


tab <- cmip6 %>% subset(period == "1985_2015"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)


knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 5.3 8.4 13.8 6.6 34.1 37.2 28.3 0.3
CESM 4.8 10.8 14.8 5.4 35.8 35.4 28.5 0.3
CMCC 5.2 7.6 10.4 5.8 29.0 42.9 27.9 0.3
CMCC-ESM2 6.2 8.5 7.7 2.5 24.9 24.5 50.3 0.3
CNRM 5.8 12.1 11.3 6.2 35.4 39.2 25.2 0.2
EC-Earth3 8.3 9.8 10.7 4.4 33.2 36.7 29.9 0.2
FGOALS 5.8 11.4 10.4 4.9 32.5 25.3 42.0 0.2
GFDL-ESM4 5.9 9.1 9.1 3.9 28.0 36.5 35.3 0.2
INM 4.0 7.7 11.3 5.2 28.2 36.5 35.0 0.3
INM-CM5 3.5 8.3 12.0 5.4 29.2 40.6 29.9 0.3
MPI 8.7 11.1 10.3 4.4 34.5 34.1 31.0 0.2
MRI 6.2 10.7 9.1 3.8 29.8 38.4 31.6 0.3
NorESM-2-MM 4.7 9.9 15.9 5.6 36.1 30.3 33.2 0.3
mean 5.7 9.6 11.3 4.9 31.5 35.2 32.9 0.3
sd 1.5 1.5 2.3 1.1 6.4 5.5 6.7 0.1

3.4 SSP 2-4.5

3.4.1 2030-2060

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP245" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")


sdmm <- cmip6 %>% subset(period == "2030_2060" & model == "SSP245" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2030_2060" & model == "SSP245"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 5.0 9.2 14.6 6.5 35.3 38.1 26.2 0.3
CESM 4.7 12.4 15.3 5.6 38.0 35.6 26.2 0.3
CMCC 4.4 8.4 11.0 6.1 29.9 44.3 25.5 0.3
CMCC-ESM2 6.1 9.1 8.0 2.8 26.0 25.3 48.4 0.3
CNRM 6.0 12.7 12.1 6.8 37.6 39.8 22.3 0.2
EC-Earth3 7.5 9.2 11.1 4.5 32.3 36.3 31.2 0.2
FGOALS 6.9 13.3 18.5 7.8 46.5 30.3 22.9 0.2
GFDL-ESM4 5.4 10.1 9.5 4.4 29.4 36.9 33.4 0.2
INM 4.6 8.4 13.8 5.9 32.7 40.4 26.7 0.3
INM-CM5 3.5 8.6 12.2 5.4 29.7 43.1 26.9 0.3
MPI 9.0 11.9 11.1 4.5 36.5 35.1 28.2 0.2
MRI 6.2 11.2 9.3 4.1 30.8 40.0 28.9 0.3
NorESM-2-MM 5.3 11.8 16.9 5.6 39.6 33.5 26.7 0.3
mean 4.0 7.3 8.7 3.7 23.7 25.6 48.0 2.6
sd 1.0 1.2 2.2 0.9 5.3 3.6 4.6 0.0

3.4.2 2070-2100

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP245" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2070_2100" & model == "SSP245" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2070_2100" & model == "SSP245"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 4.8 10.5 15.0 7.3 37.6 38.2 24.1 0.3
CESM 6.1 11.6 15.4 5.4 38.5 36.9 24.3 0.3
CMCC 4.3 8.0 12.3 6.4 31.0 46.6 22.1 0.3
CMCC-ESM2 6.2 9.5 8.5 2.8 27.0 26.8 46.0 0.3
CNRM 6.3 12.9 12.4 6.8 38.4 41.0 20.4 0.2
EC-Earth3 7.6 9.1 11.0 4.4 32.1 37.3 30.3 0.2
FGOALS 6.9 14.0 18.1 8.1 47.1 30.3 22.3 0.2
GFDL-ESM4 5.9 9.7 9.8 4.4 29.8 38.4 31.6 0.2
INM 4.8 8.2 14.1 5.8 32.9 40.6 26.2 0.3
INM-CM5 3.7 9.5 11.2 5.9 30.3 43.5 25.9 0.3
MPI 9.0 11.9 11.1 4.5 36.5 35.1 28.2 0.2
MRI 6.5 11.0 9.6 4.0 31.1 40.8 27.9 0.3
NorESM-2-MM 5.9 11.7 17.0 5.9 40.5 34.0 25.3 0.3
mean 4.2 7.4 8.8 3.8 24.2 26.2 47.0 2.6
sd 1.0 1.3 2.1 1.0 5.4 3.7 4.5 0.0

3.5 SSP 3-7.0

3.5.1 2030-2060

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP370" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2030_2060" & model == "SSP370" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2030_2060" & model == "SSP370"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 4.9 9.6 14.3 6.8 35.6 37.4 26.8 0.3
CESM 5.2 11.4 15.6 5.4 37.6 35.7 26.5 0.3
CMCC 4.2 8.3 11.3 5.9 29.7 44.1 25.9 0.3
CMCC-ESM2 3.3 8.4 9.0 3.6 24.3 26.9 48.6 0.3
CNRM 6.2 12.9 12.9 7.1 39.1 39.6 21.0 0.2
EC-Earth3 7.9 9.7 11.6 4.8 34.0 38.5 27.3 0.2
FGOALS 6.1 11.9 11.9 5.3 35.2 27.5 37.0 0.2
GFDL-ESM4 6.3 9.4 9.9 4.1 29.7 35.9 34.1 0.2
INM 4.6 8.3 13.5 5.9 32.3 40.3 27.1 0.3
INM-CM5 3.4 8.8 12.1 5.6 29.9 43.2 26.5 0.3
MPI 8.7 12.2 11.0 4.5 36.4 36.4 26.9 0.2
MRI 6.5 10.8 9.4 3.6 30.3 40.2 29.2 0.3
NorESM-2-MM 5.4 11.0 17.9 5.7 40.0 33.2 26.5 0.3
mean 3.9 7.1 8.6 3.7 23.3 25.6 48.5 2.6
sd 1.1 1.1 1.8 0.8 4.8 3.6 4.8 0.0

3.5.2 2070-2100

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP370" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2070_2100" & model == "SSP370" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2070_2100" & model == "SSP370"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 5.5 10.5 15.4 6.9 38.3 37.3 24.2 0.3
CESM 6.0 12.2 15.6 5.7 39.5 37.4 22.9 0.3
CMCC 5.2 9.1 12.5 7.3 34.1 44.1 21.6 0.3
CMCC-ESM2 6.9 9.9 8.1 3.0 27.9 25.9 45.9 0.3
CNRM 7.3 12.9 13.2 6.7 40.1 39.3 20.3 0.2
EC-Earth3 7.6 9.8 12.4 4.7 34.5 40.3 24.9 0.2
FGOALS 6.7 13.3 13.3 4.4 37.7 28.0 34.1 0.2
GFDL-ESM4 7.5 8.5 10.9 4.7 31.6 37.0 31.1 0.2
INM 6.0 8.2 14.7 6.1 35.0 39.2 25.5 0.3
INM-CM5 3.9 10.3 14.1 6.1 34.4 40.4 24.9 0.3
MPI 8.7 12.2 11.0 4.5 36.4 36.4 26.9 0.2
MRI 7.2 11.9 9.9 4.1 33.1 39.0 27.6 0.3
NorESM-2-MM 5.9 11.8 17.5 5.9 41.1 33.7 24.9 0.3
mean 4.5 7.5 9.0 3.8 24.8 25.6 47.0 2.6
sd 0.9 1.2 1.8 0.9 4.8 3.5 4.7 0.0

3.6 SSP 5-8.5

3.6.1 2030-2060

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP585" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2030_2060" & model == "SSP585" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2030_2060" & model == "SSP585"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 4.8 9.2 15.0 7.5 36.5 37.6 25.6 0.3
CESM 5.8 11.5 15.6 5.2 38.1 36.0 25.7 0.3
CMCC 4.4 8.0 11.5 6.5 30.4 44.8 24.5 0.3
CMCC-ESM2 6.5 9.3 8.0 2.5 26.3 25.4 48.0 0.3
CNRM 6.2 13.3 12.4 6.9 38.8 40.0 21.0 0.2
EC-Earth3 7.8 9.4 11.8 4.6 33.6 38.8 27.4 0.2
FGOALS 6.8 13.5 18.1 7.7 46.1 30.8 22.8 0.2
GFDL-ESM4 6.1 9.4 9.7 4.7 29.9 36.4 33.4 0.2
INM 4.6 8.4 14.0 6.0 33.0 40.6 26.2 0.3
INM-CM5 3.5 8.6 12.4 5.8 30.3 43.6 25.8 0.3
MPI 8.6 12.4 11.1 4.5 36.6 37.0 26.2 0.2
MRI 6.0 11.1 9.6 4.0 30.7 40.8 28.1 0.3
NorESM-2-MM 5.5 11.4 17.6 5.9 40.4 33.2 26.2 0.3
mean 4.1 7.3 8.9 3.8 24.1 26.0 47.3 2.6
sd 1.0 1.3 2.1 1.0 5.4 3.7 4.7 0.0

3.6.2 2070-2100

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP585" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2070_2100" & model == "SSP585" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2070_2100" & model == "SSP585"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(13,14), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 6.1 10.8 15.9 6.5 39.3 37.5 22.9 0.3
CMCC 4.4 9.3 13.5 6.7 33.9 46.7 19.2 0.3
CMCC-ESM2 5.6 10.5 8.9 2.9 27.9 28.1 43.8 0.3
CNRM 6.0 13.6 12.8 6.4 38.8 40.9 20.1 0.2
EC-Earth3 7.7 9.6 13.0 4.9 35.2 40.5 23.9 0.2
FGOALS 15.6 15.0 18.3 6.8 55.7 24.4 19.8 0.2
GFDL-ESM4 6.5 10.0 10.8 4.8 32.1 37.6 30.1 0.2
INM 4.4 8.7 14.8 6.0 33.9 41.7 24.1 0.3
INM-CM5 3.9 9.8 12.0 6.0 31.7 44.1 23.9 0.3
MPI 8.7 12.3 11.1 4.7 36.8 36.8 26.2 0.2
MRI 7.1 11.5 10.6 4.1 33.3 40.1 26.4 0.3
NorESM-2-MM 6.4 12.8 17.4 6.5 43.1 34.3 22.3 0.3
mean 4.8 7.8 9.2 3.8 25.6 26.2 45.5 2.6
sd 2.1 1.3 2.0 0.8 6.2 4.4 4.6 0.0

4 Multimodel average

Here we compute the multimodel mean and median aridity index. We compare the category obtained by using the mean aridity index (cat.mean) to the category obtained in most models (cat.maj).

4.1 Summary dataframe

4.1.1 Difference with AI ref for all SSP

Compute the difference between future AI, precipitation and temperature for future scenarios, for all models.

cmip6 <- read.table("cmip6_AI.txt")

cmip6.ref <- cmip6 %>% filter(period == "1970_2000") %>% select(c("lon","lat","model","period","source","AI","spr","t","sET0")) %>% setNames(c("lon","lat","model","period","source","AI.ref","spr.ref","t.ref","ET0.ref"))

cmip6ext <- merge(cmip6, select(cmip6.ref, c("lon","lat","source","AI.ref","spr.ref","t.ref","ET0.ref")) , by = c("lon","lat","source"), all = T) %>%  mutate(diff.AI = AI - AI.ref, diff.pr = spr - spr.ref, diff.t = t - t.ref, diff.ET0 = sET0 - ET0.ref)

write.table(cmip6ext, "cmip6ext.txt")

Create a summary table with multimodel average.

cmip6 <- read.table("cmip6ext.txt")
cmip6s <- cmip6 %>% 
  group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, period, z) %>% 
  dplyr::summarise(t.mean = mean(t, na.rm = T), t.sd = sd(t, na.rm = T), # t in °C
                   pr.mean = mean(spr, na.rm = T) , pr.sd = sd(spr, na.rm = T),  # pr in mm/y
                   sfcWind.mean = mean(sfcWind, na.rm = T), sfcWind.sd = sd(sfcWind, na.rm = T),
                   Rn.mean = mean(Rn, na.rm = T), Rn.sd = sd(Rn, na.rm = T), # in MJ/m2/day
                   ea.mean = mean(ea, na.rm = T), ea.sd = sd(ea, na.rm = T), # in kPa
                   ET0 = mean(sET0, na.rm = T), ET0.sd = sd(sET0, na.rm = T), # in mm/y
                   ET0.ref = mean(ET0.ref, na.rm = T), ET0.ref.sd = sd(ET0.ref, na.rm = T), # in mm/y
                   diff.ET0.mean = mean(diff.ET0, na.rm = T), diff.ET0.sd = sd(diff.ET0, na.rm = T),
                   AI.mean = mean(AI, na.rm = T), AI.sd = sd(AI, na.rm = T),
                   AI.med = median(AI, na.rm = T),
                   AI.ref.mean = mean(AI.ref, na.rm = T), AI.ref.sd = sd(AI.ref, na.rm = T),
                   AI.ref.med = median(AI.ref, na.rm = T),
                   diff.AI.mean = mean(diff.AI, na.rm = T), diff.AI.sd = sd(diff.AI, na.rm = T),
                   diff.AI.med = median(diff.AI, na.rm = T), 
                   spr.ref.mean = mean(spr.ref, na.rm = T), spr.ref.sd = sd(spr.ref, na.rm = T),
                   diff.pr.mean = mean(diff.pr, na.rm = T), diff.pr.sd = sd(diff.pr, na.rm = T),
                   t.ref.mean = mean(t.ref, na.rm = T), t.ref.sd = sd(t.ref, na.rm =T),
                   diff.t.mean = mean(diff.t, na.rm = T), diff.t.sd = sd(diff.t, na.rm = T),
                   Hyperarid = sum(cat.AI == "Hyperarid"),
                   Arid = sum(cat.AI == "Arid"),
                   'Semi-arid' = sum(cat.AI == "Semi-arid"),
                   'Dry subhumid' = sum(cat.AI == "Dry subhumid") ,
                   'Humid' = sum(cat.AI == "Humid"),
                   'Cold' = sum(cat.AI == "Cold"),
                   'NA' = sum(is.na(cat.AI)),
                   nb.models = n(),
                    plus.AI = table(diff.AI > 0)["TRUE"],
                    minus.AI = table(diff.AI < 0)["TRUE"]) %>%
ungroup() %>% as.data.frame()

cmip6s$cat.maj <- colnames(select(cmip6s, c("Hyperarid","Arid",'Semi-arid',"Dry subhumid",'Humid',"Cold","NA")))[max.col(select(cmip6s, c("Hyperarid","Arid",'Semi-arid',"Dry subhumid",'Humid',"Cold","NA")),ties.method = "first")]
cmip6s$nb.maj <- cmip6s %>% select(c("Hyperarid","Arid",'Semi-arid',"Dry subhumid",'Humid',"Cold","NA")) %>% apply(1, max)
cmip6s$prop.maj <- with(cmip6s, nb.maj / nb.models)



breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")

cmip6s$cat.AI <- cmip6s$AI.mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6s$cat.AI[which(cmip6s$ET0 < 400)] <- "Cold"

cmip6s$cat.AI.ref <- cmip6s$AI.ref.mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6s$cat.AI.ref[which(cmip6s$ET0.ref < 400)] <- "Cold"

cmip6s$same.catAI <- with(cmip6s, ifelse(cat.maj == cat.AI, "same", "diff"))

cmip6s$cat.dryer <- with(cmip6s, 
          ifelse(cat.AI.ref == "Arid" & cat.AI %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
          ifelse(cat.AI.ref == "Arid" & cat.AI == "Hyperarid", "dryer",
          ifelse(cat.AI.ref == "Semi-arid" & cat.AI %in% c("Cold","Humid","Dry subhumid"), "wetter",
          ifelse(cat.AI.ref == "Semi-arid" & cat.AI %in% c("Hyperarid","Arid"), "dryer",
          ifelse(cat.AI.ref == "Dry subhumid" & cat.AI %in% c("Cold","Humid"), "wetter",
          ifelse(cat.AI.ref == "Dry subhumid" & cat.AI %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.ref == "Humid" & cat.AI == "Cold","wetter",
          ifelse(cat.AI.ref == "Humid" & cat.AI %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.ref == "Cold" & cat.AI %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.ref == "Hyperarid" & cat.AI %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
          ifelse(cat.AI.ref == cat.AI, "same",
          NA))))))))))))

cmip6sref <- cmip6s %>% filter(period == "1970_2000") %>% select(c("lon","lat","model","period","cat.maj")) %>% setNames(c("lon","lat","model","period","cat.maj.ref"))

cmip6s <- merge(cmip6s, select(cmip6sref, c("lon","lat","cat.maj.ref")), 
                 by = c("lon","lat"), all.x = T) 


cmip6s$cat.maj.dryer <- with(cmip6s, 
        ifelse(cat.maj.ref == "Arid" & cat.maj %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
        ifelse(cat.maj.ref == "Arid" & cat.maj == "Hyperarid", "dryer",
        ifelse(cat.maj.ref == "Semi-arid" & cat.maj %in% c("Cold","Humid","Dry subhumid"), "wetter",
        ifelse(cat.maj.ref == "Semi-arid" & cat.maj %in% c("Hyperarid","Arid"), "dryer",
        ifelse(cat.maj.ref == "Dry subhumid" & cat.maj %in% c("Cold","Humid"), "wetter",
        ifelse(cat.maj.ref == "Dry subhumid" & cat.maj %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
        ifelse(cat.maj.ref == "Humid" & cat.maj == "Cold","wetter",
        ifelse(cat.maj.ref == "Humid" & cat.maj %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
        ifelse(cat.maj.ref == "Cold" & cat.maj %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
        ifelse(cat.maj.ref == "Hyperarid" & cat.maj %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
        ifelse(cat.maj.ref == cat.maj, "same",
        NA))))))))))))


cmip6s$Continent[which(cmip6s$Continent == "EUROPE-AFRICA")] <- "MED-BASIN"

write.table(cmip6s, "cmip6s.txt")
cmip6s <- read.table("cmip6s.txt")

4.2 Compare cat mean and cat maj

cat mean = category of aridity obtained by using the multimodel mean AI cat maj = category of aridity obtained by using the category majoritarily obtained in the 13 models

4.2.1 Where are cat mean and cat maj different


cmip6s %>% group_by(model, period) %>% summarise(maj.equal.mean = table(same.catAI)[2], maj.diff.mean = table(same.catAI)[1], prop.same = round(maj.equal.mean/n()*100,1)) %>% knitr::kable(col.names = c("Model","Period","Number of gridcells in which cat.maj = cat.mean", "Number of gridcells in which cat.maj /= cat.mean", "Proportion of cells with the same category (%)")) %>% kableExtra::kable_styling(bootstrap_options = "bordered")
Model Period Number of gridcells in which cat.maj = cat.mean Number of gridcells in which cat.maj /= cat.mean Proportion of cells with the same category (%)
SSP245 2030_2060 19884 1798 89.3
SSP245 2070_2100 19890 1792 89.3
SSP370 2030_2060 19988 1694 89.8
SSP370 2070_2100 19911 1771 89.4
SSP585 2030_2060 19799 1883 88.9
SSP585 2070_2100 19796 1886 88.9
historical 1850_1880 19907 1775 89.4
historical 1970_2000 19880 1802 89.3
historical 1985_2015 19936 1746 89.6

cmip6s$cat.AI <- factor(cmip6s$cat.AI, levels =  c("Hyperarid","Arid","Semi-arid","Dry subhumid", "Humid","Cold"))
cmip6s$cat.maj <- factor(cmip6s$cat.maj, levels =  c("Hyperarid","Arid","Semi-arid","Dry subhumid", "Humid","Cold"))

ggplot(subset(cmip6s, cat.AI != "Cold"))+
  geom_col(aes(x= period, y = cat.AI, fill = cat.AI), just = -0.2, width = 0.3)+
  geom_col(aes(x = period, y = cat.maj, fill = cat.maj), just = 1.2, width = 0.3)+
  scale_fill_manual(values = col.cat, na.translate = F)+
  labs(fill = "", y ="")+
  facet_grid(rows = vars(model), switch = "y")+
  theme_minimal()+theme(axis.text.y = element_blank())




tab <- rbind(table(cmip6s$cat.AI)/length(cmip6s$cat.AI)*100, table(cmip6s$cat.maj)/length(cmip6s$cat.maj)*100) %>% t() %>% as.data.frame() %>% setNames(c("Cat mean", "Cat maj"))
knitr::kable(tab, digits =2, caption = "Proportion of gridcells in each aridity category") %>% kable_styling(bootstrap_options = "bordered")
Proportion of gridcells in each aridity category
Cat mean Cat maj
Hyperarid 3.51 4.18
Arid 6.25 7.43
Semi-arid 7.63 9.19
Dry subhumid 3.39 0.84
Humid 28.75 28.53
Cold 47.90 47.23
map_list <- list()

for(i in unique(cmip6s$model)){
  df <- subset(cmip6s, model == i)
  for(j in unique(df$period)){
    index <- paste(i, j, sep = "")
    g <- ggplot()+geom_raster(data = subset(df, period == j & same.catAI == "diff"), aes(x = lon, y = lat, fill = cat.AI))+
      borders(colour = "grey60")+
      scale_fill_manual(values =  col.cat, na.translate = F)+
      labs(title = paste(i, j , sep = " "), fill = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }
}

col1 <- ggarrange(plotlist = map_list, nrow = 9, ncol = 1, common.legend = T, legend = "bottom") %>% annotate_figure(top = "cat.mean for gridcells in which\ncat.maj and cat.mean are different")

map_list2 <- list()

for(i in unique(cmip6s$model)){
  df <- subset(cmip6s, model == i)
  for(j in unique(df$period)){
    index <- paste(i, j, sep = "")
    g <- ggplot()+geom_raster(data = subset(df, period == j & same.catAI == "diff"), aes(x = lon, y = lat, fill = cat.maj))+
      borders(colour = "grey60")+
      scale_fill_manual(values =  col.cat, na.translate = F)+
      labs(title = paste(i, j , sep = " "), fill = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list2[[index]] <- g
  }
}

col2 <- ggarrange(plotlist = map_list2, nrow = 9, ncol = 1, common.legend = T, legend = "bottom") %>% annotate_figure(top = "cat.maj for gridcells in which\ncat.maj and cat.mean are different")

ggarrange(col1,col2,ncol = 2, common.legend = T)

4.2.2 Plot cat maj and cat mean for reference period

ggpubr::ggarrange(plotlist = list(
  
  ggplot() + 
    geom_raster(data = subset(cmip6s, period == "1970_2000"), aes(x=lon, y = lat,  fill = cat.maj))+
    geom_point(data = subset(cmip6s, period == "1970_2000" & nb.maj > 10), aes(x =lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
    scale_fill_manual(values = col.cat)+
    labs(title = "Category chosen in majority by models", fill = "")+
    ylim(-55,90)+
    theme_void(base_size = 15)+
    theme(legend.position = "none"),
  
  ggplot() + 
    geom_raster(data = subset(cmip6s, period == "1970_2000"), aes(x=lon, y = lat,  fill = cat.AI))+
   geom_point(data = subset(cmip6s, period == "1970_2000" & nb.maj > 10), aes(x =lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
    scale_fill_manual(values = col.cat)+
    labs(title = "Multimodel average categories of AI", fill = "")+
    ylim(-55,90)+
    theme_void(base_size = 15)+
    theme(legend.position = "none")),
  
  ncol = 2)
map_list <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = cat.maj))+
    geom_point(data = subset(cmip6s, period == i & model == "historical" & nb.maj > 10), aes(x= lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list[[i]] <- g
}


map_list2 <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI))+
   # geom_point(data = subset(cmip6s, period == i & model == "historical" & nb.maj > 10), aes(x= lon, y = lat), size = 0.5, col = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list2[[i]] <- g
}


ggarrange(
  annotate_figure(ggarrange(plotlist = map_list, nrow = 3, common.legend = T, legend = "bottom"), top = "Majority category among models"),
  annotate_figure(ggarrange(plotlist = map_list2, nrow = 3, common.legend = T, legend = "bottom"), top = "Multimodel average category"),
  ncol = 2, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = cat.maj))+
         geom_point(data = subset(cmip6s, period == j & model == i & nb.maj > 10), aes(x= lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
      scale_fill_manual(values = col.cat, na.translate = F)+
      labs(fill = "Cat AI maj", title = paste(i, j, sep = ", "))+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}


map_list2 <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = cat.AI))+
    #     geom_point(data = subset(cmip6s, period == j & model == i & nb.maj > 10), aes(x= lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
      scale_fill_manual(values = col.cat, na.translate = F)+
      labs(fill = "Cat AI mean", title = paste(i, j, sep = ", "))+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list2[[index]] <- g
  }}

ggarrange(plotlist = map_list, nrow = 3, ncol = 2, common.legend = T, legend = "bottom")

ggarrange(plotlist = map_list2, nrow = 3, ncol = 2, common.legend = T, legend = "bottom")

4.3 Multimodel maps of cat mean


cmip6s$cat.AI <- factor(cmip6s$cat.AI, levels = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))

map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + 
    geom_tile(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI))+
    geom_point(data = subset(cmip6s,nb.maj > 10), aes(x = lon, y = lat), shape = ".", col = "grey20", size = 0.1)+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = paste("Category AI",i , sep = " "), fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + 
      geom_tile(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = cat.AI))+
      geom_point(data = subset(cmip6s,nb.maj > 10), aes(x = lon, y = lat), shape = ".", col = "grey20", size = 0.1)+
      borders(colour = "grey60")+
      scale_fill_manual(values =  col.cat, na.translate = F)+
      labs(title = paste("Category AI",index, sep = " "), fill = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")


cmip6s$cat.AI <- factor(cmip6s$cat.AI, levels = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))

ggplot() + 
    geom_tile(data = subset(cmip6s, period == "1970_2000" & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "Category AI 1970-2000", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")


ggplot() + 
    geom_tile(data = subset(cmip6s, period == "2070_2100" & model == "SSP585"), aes(x=lon, y = lat,  fill = cat.AI))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "Category AI 2070-2100, SSP 585", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")

  map_list[[i]] <- g

4.4 Table of multimodel average category in percent

# cmip6s$cat.AI <- cmip6s$AI.mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
# cmip6s$cat.AI[which(cmip6s$ET0 < 400)] <- "Cold"
cmip6s$cat.AI.med <- cmip6s$AI.med %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6s$cat.AI.med[which(cmip6s$ET0 < 400)] <- "Cold"

table.percent <- cmip6s %>%  subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  select(period, model, cat.AI) %>% group_by(period, model, cat.AI) %>% summarise(count = n()) %>% 
  ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100,1))

table.percent.med <- cmip6s %>%  subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  select(period, model, cat.AI.med) %>% setNames(c("period","model","cat.AI")) %>% group_by(period, model, cat.AI) %>% summarise(count.med = n()) %>% 
  ungroup() %>% group_by(period, model) %>% mutate(percent.med = round(count.med/sum(count.med)*100,1))
  
table <- merge(table.percent, table.percent.med, by = c("period","model","cat.AI"), all.x = T) 
kable(table)
period model cat.AI count percent count.med percent.med
1850_1880 historical Hyperarid 746 4.8 906 5.9
1850_1880 historical Arid 1340 8.7 1450 9.4
1850_1880 historical Semi-arid 1577 10.2 1742 11.3
1850_1880 historical Dry subhumid 729 4.7 767 5.0
1850_1880 historical Humid 6385 41.3 5912 38.3
1850_1880 historical Cold 4631 30.0 NA NA
1850_1880 historical NA 36 0.2 4667 30.2
1970_2000 historical Hyperarid 739 4.8 897 5.8
1970_2000 historical Arid 1340 8.7 1450 9.4
1970_2000 historical Semi-arid 1559 10.1 1681 10.9
1970_2000 historical Dry subhumid 729 4.7 767 5.0
1970_2000 historical Humid 6174 40.0 5746 37.2
1970_2000 historical Cold 4867 31.5 NA NA
1970_2000 historical NA 36 0.2 4903 31.7
1985_2015 historical Hyperarid 732 4.7 879 5.7
1985_2015 historical Arid 1345 8.7 1431 9.3
1985_2015 historical Semi-arid 1541 10.0 1703 11.0
1985_2015 historical Dry subhumid 723 4.7 733 4.7
1985_2015 historical Humid 6062 39.3 5657 36.6
1985_2015 historical Cold 5005 32.4 NA NA
1985_2015 historical NA 36 0.2 5041 32.6
2030_2060 SSP245 Hyperarid 772 5.0 882 5.7
2030_2060 SSP245 Arid 1401 9.1 1606 10.4
2030_2060 SSP245 Semi-arid 1731 11.2 1833 11.9
2030_2060 SSP245 Dry subhumid 729 4.7 833 5.4
2030_2060 SSP245 Humid 6397 41.4 5876 38.0
2030_2060 SSP245 Cold 4378 28.3 NA NA
2030_2060 SSP245 NA 36 0.2 4414 28.6
2030_2060 SSP370 Hyperarid 717 4.6 886 5.7
2030_2060 SSP370 Arid 1483 9.6 1532 9.9
2030_2060 SSP370 Semi-arid 1686 10.9 1810 11.7
2030_2060 SSP370 Dry subhumid 676 4.4 843 5.5
2030_2060 SSP370 Humid 6394 41.4 5885 38.1
2030_2060 SSP370 Cold 4452 28.8 NA NA
2030_2060 SSP370 NA 36 0.2 4488 29.1
2030_2060 SSP585 Hyperarid 792 5.1 909 5.9
2030_2060 SSP585 Arid 1369 8.9 1569 10.2
2030_2060 SSP585 Semi-arid 1759 11.4 1903 12.3
2030_2060 SSP585 Dry subhumid 782 5.1 899 5.8
2030_2060 SSP585 Humid 6502 42.1 5924 38.4
2030_2060 SSP585 Cold 4204 27.2 NA NA
2030_2060 SSP585 NA 36 0.2 4240 27.5
2070_2100 SSP245 Hyperarid 792 5.1 920 6.0
2070_2100 SSP245 Arid 1431 9.3 1586 10.3
2070_2100 SSP245 Semi-arid 1742 11.3 1876 12.1
2070_2100 SSP245 Dry subhumid 777 5.0 878 5.7
2070_2100 SSP245 Humid 6532 42.3 6014 38.9
2070_2100 SSP245 Cold 4134 26.8 NA NA
2070_2100 SSP245 NA 36 0.2 4170 27.0
2070_2100 SSP370 Hyperarid 899 5.8 992 6.4
2070_2100 SSP370 Arid 1399 9.1 1594 10.3
2070_2100 SSP370 Semi-arid 1842 11.9 1935 12.5
2070_2100 SSP370 Dry subhumid 822 5.3 918 5.9
2070_2100 SSP370 Humid 6324 40.9 5847 37.9
2070_2100 SSP370 Cold 4122 26.7 NA NA
2070_2100 SSP370 NA 36 0.2 4158 26.9
2070_2100 SSP585 Hyperarid 852 5.5 964 6.2
2070_2100 SSP585 Arid 1412 9.1 1614 10.5
2070_2100 SSP585 Semi-arid 1857 12.0 1941 12.6
2070_2100 SSP585 Dry subhumid 805 5.2 905 5.9
2070_2100 SSP585 Humid 6686 43.3 6188 40.1
2070_2100 SSP585 Cold 3796 24.6 NA NA
2070_2100 SSP585 NA 36 0.2 3832 24.8

ggplot(subset(table, period == "1970_2000"), aes(x = cat.AI))+
  geom_col(aes(y = percent, fill = "mean"), position = position_nudge(x = - 0.2), width = 0.4)+
  geom_col(aes(y = percent.med, fill = "med"), position = position_nudge(x = + 0.2), width = 0.4)+
  scale_fill_manual(values = c("#ABD9E9", "#FDAE61"))+
  labs(x="Aridity category", fill = "Using mean or median")+
  theme_minimal()

5 Maps multimodel AI

quantile(cmip6s$AI.mean, probs = seq(0,1,0.1), na.rm = T) 
##           0%          10%          20%          30%          40%          50% 
## 2.902197e-03 1.993698e-01 5.922059e-01 1.019621e+00 1.457766e+00 2.127218e+00 
##          60%          70%          80%          90%         100% 
## 3.578980e+00 1.168470e+01 2.438867e+01 4.856852e+01 5.329875e+06

ai.breaks <- c(0,0.03,0.2,0.5,0.65,1,10,20,30)
colscale <- brewer.pal(n = 11, name ="RdYlBu")[-c(5,6)]
map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = AI.mean))+
    borders(colour = "grey60")+
    binned_scale(aesthetics = "fill", breaks = ai.breaks, palette = function(x) colscale,
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "Aridity index")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = AI.mean))+
      borders(colour = "grey60")+
      binned_scale(aesthetics = "fill", breaks = ai.breaks, palette = function(x) colscale,
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Aridity index", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

6 AI changes over time and scenarios

v<-NULL
v2<-NULL
for (N in 5:30) {
  v<-c(v,sum(sapply(floor(2*N/3):N,function(i) choose(N,i)*(1/2)^N*2)))# proba d'avoir 2/3 de même signe à partir de N modèles, lorsque N change
  v2<-c(v2,sum(sapply(floor(2*N/3):N,function(i) choose(N,i)*(1/2)^N)))# proba d'avoir 2/3 de signe >0 à partir de N modèles, lorsque N change
}


plot(5:30,v,type="l",ylab="Proba(>2/3 models...)",xlab="Nb of models")
lines(5:30,v2,col="red")

6.1 Changes AI for summary CMIP6

6.2 Boxplot of changes AI

ggpubr::ggarrange(plotlist = list(
  ggplot(subset(cmip6s, diff.AI.mean < 1 & diff.AI.mean > -1 & model == "historical"))+
    geom_boxplot(aes(x=period, y = - diff.AI.mean, col = period), outliers = F)+
    scale_color_manual(values = c("#D73027","#FDAE61","#FEE08B"))+
    labs(x = "", y = "Difference of AI between each period and\nthe reference period 1970-2000")+
    facet_grid(cols = vars(model))+
    theme_minimal() + theme(legend.position = "none", axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)),
  ggplot(subset(cmip6s, diff.AI.mean < 1 & diff.AI.mean > -1 & model != "historical"))+
    geom_boxplot(aes(x=period, y = diff.AI.mean, col = period), outliers = F)+
    scale_color_manual(values = c("#ABD9E9","#4575B4"))+
    scale_y_continuous(position = "right")+
    labs(y = "", x = "")+
    facet_grid(cols = vars(model))+
    theme_minimal() + theme(legend.position = "none", axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))),
  widths = c(1.5,3),
  ncol = 2, nrow = 1
)

6.3 Strictly wetter or dryer (reference period: 1970-2000)

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
  index <- paste(i, j, sep = ", ")
  g <- ggplot()+
  geom_point(data = subset(cmip6s, model == i & period == j & plus.AI > 10), aes(x = lon, y = lat, col = "plus"), shape = "+")+
  geom_point(data = subset(cmip6s, model == i & period == j & minus.AI > 10), aes(x = lon, y = lat, col = "minus"), shape = "-")+
    scale_color_manual(values = c("#ABD9E9", "#D73027"))+
    labs(color = "", x = "", y = "", title = index)+
  theme_minimal()
  map_list[[index]] <- g
  }
}

ggarrange(plotlist = map_list, ncol = 2, nrow = 3, common.legend = T)

map_list <- list()

for(i in c("1850_1880","1985_2015")){
    index <- paste(i, j, sep = " , ")
  g <- ggplot() + geom_tile(data = subset(cmip6s, period == i & model == "historical" & diff.AI.mean > 0), aes(x=lon, y = lat,  fill = "wetter"))+
    geom_tile(data = subset(cmip6s, period == i & model == "historical" & diff.AI.mean < 0), aes(x=lon, y = lat,  fill = "dryer"))+
    geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = "+", size = 0.5)+
     geom_point(data = subset(cmip6s, period == i & model == "historical" & minus.AI > 10), aes(x=lon, y = lat), shape = "-", size = 0.5)+
    borders(colour = "grey60")+
    scale_fill_manual(values =  c("#FDAE61", "#ABD9E9"), na.translate = F)+
    labs(title = paste("AI changes between",i , "and 1970-2000", sep = " "), fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list[[i]] <- g
  }

ggarrange(plotlist = map_list, ncol = 2)


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & diff.AI.mean > 0), aes(x=lon, y = lat,  fill = "wetter"))+
      geom_tile(data = subset(cmip6s, period == j & model == i & diff.AI.mean < 0), aes(x=lon, y = lat,  fill = "dryer"))+
       geom_point(data = subset(cmip6s, period == j & model == i & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s, period == j & model == i & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
      borders(colour = "grey60")+
      scale_fill_manual(values = c("#F46D43", "#74ADD1"), na.translate = F)+
      labs(fill = "Compared to 1970-2000", title = i)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}



ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

6.4 Wetter or dryer, binned (reference period: 1970-2000)

ano.quantiles <- cmip6s %>% group_by(model, period) %>% summarise(q10 = quantile(diff.AI.mean, probs = 0.1, na.rm =T), 
                                                                  q25 = quantile(diff.AI.mean, probs = 0.25, na.rm =T), 
                                                                  q50 = quantile(diff.AI.mean, probs = 0.5, na.rm =T), 
                                                                  q75 = quantile(diff.AI.mean, probs = 0.75, na.rm =T), 
                                                                  q90 = quantile(diff.AI.mean, probs = 0.9, na.rm =T)) 

ano.quantiles <- cmip6s %>% summarise(q10 = quantile(diff.AI.mean, probs = 0.1, na.rm =T), 
                                      q25 = quantile(diff.AI.mean, probs = 0.25, na.rm =T), 
                                      q50 = quantile(diff.AI.mean, probs = 0.5, na.rm =T), 
                                      q75 = quantile(diff.AI.mean, probs = 0.75, na.rm =T), 
                                      q90 = quantile(diff.AI.mean, probs = 0.9, na.rm =T)) 
qbreaks <- c(-5, -1,-0.1, -0.01, 0, 0.01, 0.1, 1,5)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")

“diff.AI.mean” is the multimodel AI average minus the reference AI (period 1970-2000). More blue: AI mean is superior to AI ref, hence wetter than the ref. More red: the anomalie is negative, hence AI mean is inferior to AI ref, hence dryer.

AS OF NOW THE SAME QUANTILES ARE USED FOR ALL MAPS.

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  #quant = subset(ano.quantiles, model == "historical" & period == i)
  #b = quant[1,c(3:7)] %>% as.numeric
  b = qbreaks
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = diff.AI.mean))+
    borders(colour = "grey60")+
           geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s,  period == i & model == "historical"  & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
    binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale),
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    b = qbreaks
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = diff.AI.mean))+
      borders(colour = "grey60")+
                 geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s,  period == i & model == "historical"  & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
      binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale),
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Compared to 1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}



ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

6.5 Wetter or dryer, binned, % (reference period: 1970-2000)

cmip6s$diff.AI.mean.percent <- with(cmip6s, (AI.mean-AI.ref.mean)/AI.ref.mean*100)


pbreaks <- c(-40, -30,-20, -10, 0, 10, 20, 30, 40)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef","white","white", "#fddbc7", "#f4a582", "#d6604d", "#b2182b")

“diff.AI.mean” is the multimodel AI average minus the reference AI (period 1970-2000). More blue: AI mean is superior to AI ref, hence wetter than the ref. More red: the anomalie is negative, hence AI mean is inferior to AI ref, hence dryer.

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  #quant = subset(ano.quantiles, model == "historical" & period == i)
  #b = quant[1,c(3:7)] %>% as.numeric
  b = pbreaks
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = diff.AI.mean.percent))+
    borders(colour = "grey60")+
           geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s,  period == i & model == "historical"  & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
    binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale))+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom") %>% annotate_figure(top = "% change of AI index compared to 1970-2000")


#,guide = guide_legend(label.theme = element_text(angle = 0))

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    b = pbreaks
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = diff.AI.mean.percent))+
      borders(colour = "grey60", size = 0.5)+
                 geom_point(data = subset(cmip6s, period == j & model == i & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s,  period == j & model == i & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
      binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale))+
      labs(fill = "% change of AI index\ncompared to 1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}


ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")


ggsave("map_AIpercent_SSP585_2100.png", map_list[[index]])

7 Maps of difference between future and reference category of AI, using cat maj

7.1 Detailed changes and colors cat maj

cmip6s$cat.maj[which(is.na(cmip6s$cat.maj) == T)] <- "NA"
cmip6s$cat.maj.ref[which(is.na(cmip6s$cat.maj.ref) == T)] <- "NA"
cmip6s$change.catmaj <- with(cmip6s, ifelse(cat.maj == cat.maj.ref, 0, 1))

cmip6s$change.cat.maj <- with(cmip6s, paste(cat.maj.ref, cat.maj, sep = " to "))

cmip6s$change.cat.maj <- factor(cmip6s$change.cat.maj, 
                                levels = c("Cold to Humid", 
                                           "Cold to Dry subhumid",
                                           "Cold to Semi-arid",
                                           "Humid to Arid",
                                           "Humid to Semi-arid",
                                           "Humid to Dry subhumid",
                                           "Dry subhumid to Semi-arid",
                                           "Dry subhumid to Arid",
                                           "Semi-arid to Arid",
                                           "Arid to Hyperarid", 
                                           "Hyperarid to Arid",
                                           "Arid to Semi-arid",
                                           "Arid to Dry subhumid",
                                           "Arid to Humid",
                                           "Semi-arid to Dry subhumid",
                                           "Semi-arid to Humid",
                                           "Semi-arid to Cold",
                                           "Dry subhumid to Humid",
                                           "Dry subhumid to Cold",
                                           "Humid to Cold"))

colors_changes <- c("Cold to Humid" = "#FEE08B", 
                                           "Cold to Dry subhumid" = "#FDCB6A",
                                           "Cold to Semi-arid" = "#FDB463",
                                           "Humid to Arid" = "#FDAE61",
                                           "Humid to Semi-arid" = "#F68D54",
                                           "Humid to Dry subhumid" =  "#F46B43",
                                           "Dry subhumid to Semi-arid" = "#F03B20",
                                           "Dry subhumid to Arid" = "#E63125",
                                           "Semi-arid to Arid" = "#D73027",
                                           "Arid to Hyperarid" = "#BD0026", 
                                           "Hyperarid to Arid" = "#ABD9E9",
                                           "Arid to Semi-arid" = "#91C4DC",
                                           "Arid to Dry subhumid" = "#7FB2D0",
                                           "Arid to Humid" = "#6A9EC4",
                                           "Semi-arid to Dry subhumid" = "#5893B8",
                                           "Semi-arid to Humid" = "#4575B4",
                                           "Semi-arid to Cold" = "#5D83C9",
                                           "Dry subhumid to Humid" = "#77A1E2",
                                           "Dry subhumid to Cold" = "#A7C7E9",
                                           "Humid to Cold" = "#DAE2F7")

7.2 Changes towards dryer categories cat maj

7.2.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + 
    geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & cat.maj.dryer == "dryer"), 
                aes(x=lon, y = lat,  fill = change.cat.maj))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 3))
  map_list[[i]] <- g
}


ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

7.2.1.1 Future


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmaj == 1 & cat.maj.dryer == "dryer"), 
                              aes(x=lon, y = lat,  fill = change.cat.maj))+
      borders(colour = "grey60")+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 3))
    map_list[[index]] <- g
  }}

ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

7.3 Changes towards wetter categories cat maj

7.3.0.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & cat.maj.dryer == "wetter"), 
                              aes(x=lon, y = lat,  fill = change.cat.maj))+
    borders(colour = "grey60", ylim = c(-40,40))+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+
    ylim(-50,50)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 3))
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

7.3.0.2 Future


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmaj == 1 & cat.maj.dryer == "wetter"), 
                              aes(x=lon, y = lat,  fill = change.cat.maj))+
      borders(colour = "grey60", ylim = c(-40,40))+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+ylim(-50,50)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 3))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list, nrow = 3, ncol = 2, common.legend = T, legend = "bottom")

8 Map of difference between future and reference category of AI, using cat mean

8.1 Detailed changes and colors

cmip6s$change.catmean <- with(cmip6s, ifelse(cat.AI == cat.AI.ref, 0, 1))

cmip6s$change.cat.mean <- with(cmip6s, paste(cat.AI.ref, cat.AI, sep = " to "))
unique(cmip6s$change.cat.mean)
##  [1] "NA to NA"                     "Cold to Cold"                
##  [3] "Semi-arid to Semi-arid"       "Semi-arid to Arid"           
##  [5] "Arid to Arid"                 "Arid to Hyperarid"           
##  [7] "Hyperarid to Hyperarid"       "Hyperarid to Arid"           
##  [9] "Dry subhumid to Dry subhumid" "Dry subhumid to Semi-arid"   
## [11] "Humid to Semi-arid"           "Dry subhumid to Humid"       
## [13] "Humid to Humid"               "Humid to Dry subhumid"       
## [15] "Semi-arid to Dry subhumid"    "Arid to Semi-arid"           
## [17] "Cold to Humid"                "Humid to Cold"

table(cmip6s$change.cat.mean)
## 
##                 Arid to Arid            Arid to Hyperarid 
##                        11080                          639 
##            Arid to Semi-arid                 Cold to Cold 
##                          340                        95818 
##                Cold to Humid Dry subhumid to Dry subhumid 
##                         4393                         4022 
##        Dry subhumid to Humid    Dry subhumid to Semi-arid 
##                          288                         2242 
##                Humid to Cold        Humid to Dry subhumid 
##                          148                         2484 
##               Humid to Humid           Humid to Semi-arid 
##                        52925                          198 
##            Hyperarid to Arid       Hyperarid to Hyperarid 
##                          240                         6402 
##                     NA to NA            Semi-arid to Arid 
##                         5130                         1209 
##    Semi-arid to Dry subhumid       Semi-arid to Semi-arid 
##                          277                        12514
cmip6s$change.cat.mean <- factor(cmip6s$change.cat.mean, 
                                  levels = c("Hyperarid to Arid",
                                             "Arid to Semi-arid",
                                             "Semi-arid to Dry subhumid",
                                             "Dry subhumid to Humid",
                                             "Humid to Cold",
                                             "Cold to Humid",
                                             "Humid to Dry subhumid",
                                             "Humid to Semi-arid",
                                             "Dry subhumid to Semi-arid",
                                             "Semi-arid to Arid",
                                             "Arid to Hyperarid"))

wtd <- colorRampPalette(c("#FEE08B","firebrick"))
col.wtd <- wtd(6)

dtw <- colorRampPalette(c( "#DAE2F7","#4575B4" ))
col.dtw <- dtw(5)

colors_changes <- c("Hyperarid to Arid" = col.dtw[1],
                    "Arid to Semi-arid" = col.dtw[2],
                    "Semi-arid to Dry subhumid" = col.dtw[3],
                    "Dry subhumid to Humid" = col.dtw[4],
                    "Humid to Cold" = col.dtw[5],
                    "Cold to Humid" = col.wtd[1],
                    "Humid to Dry subhumid" = col.wtd[2],
                    "Humid to Semi-arid" = col.wtd[3],
                    "Dry subhumid to Semi-arid" = col.wtd[4],
                    "Semi-arid to Arid" = col.wtd[5],
                    "Arid to Hyperarid" = col.wtd[6])
   

8.2 Changes towards dryer categories cat mean

8.2.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + 
    geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmean == 1 & cat.dryer == "dryer"), 
                aes(x=lon, y = lat,  fill = change.cat.mean))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

8.2.2 Future


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmean == 1 & cat.dryer == "dryer"), 
                              aes(x=lon, y = lat,  fill = change.cat.mean))+
      borders(colour = "grey60")+
      scale_fill_manual(values = colors_changes)+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggarrange(plotlist = map_list, ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

8.3 Changes towards wetter categories cat mean

8.3.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + geom_tile(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & cat.dryer == "wetter"), 
                            aes(x=lon, y = lat,  fill = change.cat.mean))+
    borders(colour = "grey60", ylim = c(-40,40))+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+ylim(-50,50)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

8.3.2 Future


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmean == 1 & cat.dryer == "wetter"), 
                              aes(x=lon, y = lat,  fill = change.cat.mean))+
      borders(colour = "grey60", ylim = c(-40,40))+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

8.4 Big map of SSP585 2070-2100


g <- ggplot() + 
  geom_tile(data = subset(cmip6s, period == "2070_2100" & model == "SSP585" & change.catmean == 1), aes(x=lon, y = lat,  fill = change.cat.mean))+
      borders(colour = "grey60")+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      labs(fill = "Change of aridity category\nfor the period 2070-2100\nand SSP 5-8.5,\ncompared to 1970-2000", title = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 4))
   
print(g)

9 Tables

9.1 Changes of aridity category using cat mean

tab <- cmip6s %>% group_by(model, period, Continent) %>% summarise(Cold = table(cat.AI)[1], Hyperarid = table(cat.AI)[2], Arid = table(cat.AI)[3], Semiarid = table(cat.AI)[4], Drysubhumid = table(cat.AI)[5], Humid = table(cat.AI)[6], total = n()) 
write.table(tab, "tab.catAI.txt")

tab.future <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(period, model, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(period + model ~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

kable(tab.future[order(tab.future$model, decreasing = F),]) %>% kable_styling(bootstrap_options = "bordered") %>%
  column_spec(c(7,10), italic = T, include_thead = T) %>% row_spec(2, bold = T) 
model period Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
1 historical 1850_1880 3.4 6.0 7.1 3.3 19.8 28.8 48.9 2.6
2 historical 1970_2000 3.3 6.0 7.0 3.3 19.6 27.8 49.9 2.6
3 historical 1985_2015 3.3 6.1 6.9 3.3 19.6 27.3 50.6 2.6
4 SSP245 2030_2060 3.5 6.3 7.8 3.3 20.9 28.8 47.7 2.6
7 SSP245 2070_2100 3.6 6.4 7.8 3.5 21.3 29.4 46.6 2.6
5 SSP370 2030_2060 3.2 6.7 7.6 3.0 20.5 28.8 48.1 2.6
8 SSP370 2070_2100 4.0 6.3 8.3 3.7 22.3 28.5 46.6 2.6
6 SSP585 2030_2060 3.6 6.2 7.9 3.5 21.2 29.3 47.0 2.6
9 SSP585 2070_2100 3.8 6.4 8.4 3.6 22.2 30.1 45.1 2.6

# ggplot(tab.future, aes(x = period))+
#   geom_point(aes(y = Hyperarid, col = "Hyperarid", shape = model))+
#   geom_line(aes(y=Hyperarid, col = "Hyperarid", group = model, lty = model))+
#   geom_point(aes(y = Arid, col = "Arid", shape = model))+
#   geom_line(aes(y=Arid, col = "Arid", group = model, lty = model))+
#   geom_point(aes(y = get('Semi-arid'), col = "Semi-arid", shape = model))+
#   geom_line(aes(y=get('Semi-arid'), col = "Semi-arid", group = model, lty = model))+
#   geom_point(aes(y = get('Dry subhumid'), col = "Dry subhumid", shape = model))+
#   geom_line(aes(y=get('Dry subhumid'), col = "Dry subhumid", group = model, lty = model))+
#   geom_point(aes(y = Humid, col = "Humid", shape = model))+
#   geom_line(aes(y= Humid, col = "Humid", group = model, lty = model))+
#   scale_color_manual(values = col.cat)+
#   theme_minimal()

cmip6$cat.AI <- factor(cmip6$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))
tab.percent <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() %>% 
  group_by(period, model, cat.AI) %>% 
  summarise(mmmean = mean(count, na.rm = T), mmsd = sd(count, na.rm = T)) %>%
  mutate(percent = round(mmmean/sum(mmmean)*100, 1), sd.percent = round(mmsd/sum(mmmean)*100,1)) 
write.table(tab.percent, "tab.percent.txt")
tab.percent <- read.table("tab.percent.txt")

df245 <- rbind(subset(tab.percent, model == "historical" & period %in% c("1970_2000")), subset(tab.percent, model == "SSP245")) %>%
  group_by(period, model) %>% 
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

df245$cat.AI <- factor(df245$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

b245 <- ggplot(data = df245, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = paste(percent, " ± ", sd.percent, sep = "")), size = 3)+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 2-4.5", y = "%", x = "", col = "", fill = "")+
  theme_minimal()

df370 <- rbind(subset(tab.percent, model == "historical" & period %in% c( "1970_2000")), subset(tab.percent, model == "SSP370")) %>%
  group_by(period, model) %>% 
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

df370$cat.AI <- factor(df370$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

b370 <- ggplot(data = df370, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = paste(percent, " ± ", sd.percent, sep = "")), size = 3)+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 3-7.0", y = "%", x = "", col = "", fill = "")+
  theme_minimal()


df585 <- rbind(subset(tab.percent, model == "historical" & period %in% c("1970_2000")), subset(tab.percent, model == "SSP585")) %>%
  group_by(period, model) %>% 
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

df585$cat.AI <- factor(df585$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

b585 <- ggplot(data = df585, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = paste(percent, " ± ", sd.percent, sep = "")), size = 3)+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 5-8.5", y = "%", x = "", col = "", fill = "")+
  theme_minimal()


ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)

9.2 Changes of aridity category using cat maj

tab.maj <- cmip6s %>% group_by(model, period, Continent) %>% summarise(Cold = table(cat.maj)[1], Hyperarid = table(cat.maj)[2], Arid = table(cat.maj)[3], Semiarid = table(cat.maj)[4], Drysubhumid = table(cat.maj)[5], Humid = table(cat.maj)[6], total = n()) 
write.table(tab.maj, "tab.catmaj.txt")

tab.future.maj <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(period, model, cat.maj) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(period + model ~cat.maj) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

kable(tab.future.maj[order(tab.future.maj$model, decreasing = F),]) %>% kable_styling(bootstrap_options = "bordered") %>%
  column_spec(c(7,10), italic = T, include_thead = T) %>% row_spec(2, bold = T) 
model period Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
1 historical 1850_1880 4.2 7.1 8.7 0.8 20.8 28.3 48.3 2.6
2 historical 1970_2000 4.1 7.1 8.6 0.7 20.5 28.1 48.9 2.6
3 historical 1985_2015 4.0 7.0 8.4 0.9 20.3 27.3 49.9 2.6
4 SSP245 2030_2060 4.0 7.8 9.3 0.8 21.9 28.4 47.2 2.6
7 SSP245 2070_2100 4.2 7.7 9.4 0.9 22.2 29.1 46.1 2.6
5 SSP370 2030_2060 4.1 7.2 9.4 0.8 21.5 28.6 47.3 2.6
8 SSP370 2070_2100 4.5 7.7 10.0 0.9 23.1 28.4 45.9 2.6
6 SSP585 2030_2060 4.2 7.5 9.6 1.0 22.3 28.8 46.3 2.6
9 SSP585 2070_2100 4.5 7.9 9.4 0.9 22.7 30.0 44.6 2.6
tab.flow.catmaj <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, cat.maj) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL)

df245 <- rbind(subset(tab.flow.catmaj, model == "historical" & period %in% c("1970_2000")), subset(tab.flow.catmaj, model == "SSP245")) %>%
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b245 <- ggplot(data = df245, aes(x = period, y = percent))+
  geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 2-4.5", y = "%", x = "")+
  theme_minimal()

df370 <- rbind(subset(tab.flow.catmaj, model == "historical" & period %in% c("1970_2000")), subset(tab.flow.catmaj, model == "SSP370")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b370 <- ggplot(data = df370, aes(x = period, y = percent))+
  geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 3-7.0", y = "%", x = "")+
  theme_minimal()


df585 <- rbind(subset(tab.flow.catmaj, model == "historical" & period %in% c("1970_2000")), subset(tab.flow.catmaj, model == "SSP585")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b585 <- ggplot(data = df585, aes(x = period, y = percent))+
  geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 5-8.5", y = "%", x = "")+
  theme_minimal()


ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)

9.3 Changes by continent, cat mean


cmip6$cat.AI <- factor(cmip6$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

tab.percent.cont <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, Continent, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() %>% 
  group_by(period, model, Continent, cat.AI) %>% 
  summarise(mmmean = mean(count, na.rm = T), mmsd = sd(count, na.rm = T)) %>%
  mutate(percent = round(mmmean/sum(mmmean)*100, 1), sd.percent = round(mmsd/sum(mmmean)*100,1)) 
write.table(tab.percent.cont, "tab.percent.continent.txt")
tab.percent.cont <- read.table("tab.percent.continent.txt")

tab_list_percent <- list()
tab_list_sd <- list()

for(i in unique(tab.percent.cont$Continent)){
  tab.i <- subset(tab.percent.cont, Continent == i)
  
  df.percent <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "percent")
  
  df.sd <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "sd.percent")
  # sd for a sum: square root of sum of squared sd
  
  drycats <- which(names(df.percent) %in% c("Hyperarid","Arid","Semi-arid","Dry subhumid"))
  
  df.percent <- df.percent %>% mutate("Sum drylands" = rowSums(.[drycats], na.rm = T))
  df.sd <- df.sd %>% mutate("Sum drylands" = round(sqrt(rowSums((.[drycats])^2, na.rm = T)),2))
  
  missing <- setdiff(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"), names(df.percent))
  
  df.percent[, missing] <- 0
  df.sd[,missing] <- 0
  
  df.percent <- df.percent %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.percent[is.na(df.percent)] <- 0
  
  tab_list_percent[[i]] <- df.percent
  
  df.sd <- df.sd %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.sd[is.na(df.sd)] <- 0
  
  tab_list_sd[[i]] <- df.sd
}

k_list <- list()


kdiff_list <- list()

for(i in names(tab_list_percent)){
  
  df <- tab_list_percent[[i]] %>% select(c("model", "period"))
  
  
  
  prop.ref <- tab_list_percent[[i]] %>% filter(period == "1970_2000")
  prop.ref <- prop.ref %>% rbind(prop.ref[rep(1,dim(df)[1]-1),])
  df2 <- tab_list_percent[[i]][,c(3:10)] - prop.ref[,c(3:10)] %>% round(1)
  diff.prop.percent <- cbind(df, df2)
  
  sd.ref <- tab_list_sd[[i]] %>% filter(period == "1970_2000")
  sd.ref <- sd.ref %>% rbind(sd.ref[rep(1,dim(df)[1]-1),])
  df3 <- sqrt(tab_list_sd[[i]][,c(3:10)]^2 + sd.ref[,c(3:10)]^2) %>% round(1)
  diff.sd.percent <- cbind(df, df3)
  
  df$Hyperarid <- paste(tab_list_percent[[i]]$Hyperarid, tab_list_sd[[i]]$Hyperarid, sep = " ± ")
  df$Arid <- paste(tab_list_percent[[i]]$Arid, tab_list_sd[[i]]$Arid, sep = " ± ")
  df$'Semi-Arid' <- paste(tab_list_percent[[i]][,5], tab_list_sd[[i]][,5], sep = " ± ")
  df$'Dry subhumid' <- paste(tab_list_percent[[i]][,6], tab_list_sd[[i]][,6], sep = " ± ")
  df$'Sum drylands' <- paste(tab_list_percent[[i]][,7], tab_list_sd[[i]][,7], sep = " ± ")
  df$Humid <- paste(tab_list_percent[[i]]$Humid, tab_list_sd[[i]]$Humid, sep = " ± ")
  df$Cold <- paste(tab_list_percent[[i]]$Cold, tab_list_sd[[i]]$Cold, sep = " ± ")
  
  df.diff <- df
  df.diff$Hyperarid <- paste(round(diff.prop.percent$Hyperarid, 1), round(diff.sd.percent$Hyperarid, 1), sep = " ± ")
  df.diff$Arid <- paste(round(diff.prop.percent$Arid, 1), round(diff.sd.percent$Arid, 1), sep = " ± ")
  df.diff$'Semi-Arid' <- paste(round(diff.prop.percent[,5], 1), round(diff.sd.percent[,5], 1), sep = " ± ")
  df.diff$'Dry subhumid' <- paste(round(diff.prop.percent[,6], 1), round(diff.sd.percent[,6], 1), sep = " ± ")
  df.diff$'Sum drylands' <- paste(round(diff.prop.percent[,7],1), round(diff.sd.percent[,7], 1), sep = " ± ")
  df.diff$Humid <- paste(round(diff.prop.percent$Humid, 1), round(diff.sd.percent$Humid, 1), sep = " ± ")
  df.diff$Cold <- paste(round(diff.prop.percent$Cold, 1), round(diff.sd.percent$Cold, 1), sep = " ± ")
  
  k_list[[i]] <- kable(df[order(df$model, decreasing = F),], caption = i) %>% kable_styling(bootstrap_options = "bordered") %>%
    column_spec(8, italic = T, background = "#DAE2F7") %>% row_spec(2, bold = T) 
  kdiff_list[[i]] <- kable(df.diff[order(df.diff$model, decreasing = F),], caption = i) %>% kable_styling(bootstrap_options = "bordered") %>% 
    column_spec(8, italic = T, background = "#DAE2F7") %>% row_spec(2, bold = T) 
  
}

9.3.1 Proportion of aridity categories by continent

k_list
$AFRICA
AFRICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 27.2 ± 4.7 17.1 ± 4.9 15.6 ± 3.6 6.3 ± 1.4 66.2 ± 7.81 33.7 ± 7.1 0 ± 0
2 historical 1970_2000 26.9 ± 4.7 16.7 ± 4.2 15.1 ± 3.4 6.2 ± 1.5 64.9 ± 7.32 32.8 ± 6.2 2.3 ± 0
3 historical 1985_2015 26.5 ± 4.8 17.1 ± 4.6 15.4 ± 3.3 6 ± 1.5 65 ± 7.57 32.7 ± 6.3 2.2 ± 0
4 SSP245 2030_2060 27 ± 4.8 18 ± 4.9 16.5 ± 3.5 6.4 ± 1.7 67.9 ± 7.89 32 ± 7.6 0 ± 0
7 SSP245 2070_2100 27.9 ± 4.7 17.4 ± 5 17.2 ± 3.4 6.2 ± 1.5 68.7 ± 7.8 31.3 ± 7.9 0 ± 0
5 SSP370 2030_2060 26 ± 5 16.6 ± 3.6 17.2 ± 1.8 6.8 ± 1 66.6 ± 6.5 32.3 ± 6.7 1.1 ± 0
8 SSP370 2070_2100 28 ± 3.5 17.6 ± 4.8 16.6 ± 3.1 6.2 ± 1.6 68.4 ± 6.89 30.4 ± 6.4 1.2 ± 0
6 SSP585 2030_2060 27.1 ± 4.6 17.6 ± 4.8 16.8 ± 3.5 6.3 ± 1.5 67.8 ± 7.66 32.1 ± 7.5 0 ± 0
9 SSP585 2070_2100 29.7 ± 10.5 17.3 ± 5.3 16.2 ± 4.2 6 ± 1.5 69.2 ± 12.58 30.9 ± 8.6 0 ± 0
$ASIA
ASIA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 3.9 ± 2.3 12.2 ± 1.9 9.2 ± 1.5 4 ± 1.4 29.3 ± 3.62 32.4 ± 7.7 38 ± 9.2
2 historical 1970_2000 3.8 ± 2.1 12 ± 1.8 8.9 ± 1.5 3.8 ± 1.1 28.5 ± 3.33 30.7 ± 8.5 40.5 ± 9.6
3 historical 1985_2015 3.7 ± 2.2 11.8 ± 1.6 8.6 ± 1.5 3.7 ± 1.1 27.8 ± 3.3 30.1 ± 7.6 41.8 ± 8.8
4 SSP245 2030_2060 3.5 ± 1.8 12.5 ± 1.6 9.4 ± 1.6 4.3 ± 1.7 29.7 ± 3.35 34.3 ± 8.1 35.8 ± 9.2
7 SSP245 2070_2100 3.7 ± 2 12.5 ± 1.7 9.5 ± 1.7 4.2 ± 1.7 29.9 ± 3.56 36.2 ± 8.6 33.7 ± 9.3
5 SSP370 2030_2060 3.4 ± 2.3 12.4 ± 1.5 8.9 ± 1.7 3.7 ± 1.3 28.4 ± 3.48 34 ± 7.2 37.3 ± 9.5
8 SSP370 2070_2100 4.6 ± 1.9 12.5 ± 1.5 9.6 ± 1.5 3.8 ± 1.2 30.5 ± 3.09 35.2 ± 9.2 34.1 ± 9.7
6 SSP585 2030_2060 3.8 ± 2 12.5 ± 1.6 9.5 ± 1.8 4.2 ± 1.6 30 ± 3.52 35.4 ± 8.4 34.4 ± 9.3
9 SSP585 2070_2100 4.8 ± 3.2 13.1 ± 3.2 10.4 ± 3.4 4.2 ± 1.5 32.5 ± 5.86 36.5 ± 11.7 30.7 ± 9.9
$CENTRAL-AMERICA
CENTRAL-AMERICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0 8.4 ± 3.6 30.1 ± 5.1 16.1 ± 8.7 54.6 ± 10.71 45 ± 6.1 0 ± 0
2 historical 1970_2000 0 ± 0 8.4 ± 4.1 29.7 ± 7.6 14.7 ± 7.8 52.8 ± 11.64 44.9 ± 12 2 ± 0
3 historical 1985_2015 0 ± 0 9 ± 4.2 30.7 ± 7 14.8 ± 7.4 54.5 ± 11.02 43.8 ± 11.5 1.4 ± 0
4 SSP245 2030_2060 0 ± 0 11.9 ± 6.2 36.3 ± 8 13.9 ± 6.2 62.1 ± 11.87 37.6 ± 8.6 0 ± 0
7 SSP245 2070_2100 0 ± 0 12.3 ± 5.5 38.1 ± 8.5 14.4 ± 7.9 64.8 ± 12.84 34.8 ± 8.6 0 ± 0
5 SSP370 2030_2060 0 ± 0 12.6 ± 6.3 31.4 ± 7.1 15 ± 7.9 59 ± 12.35 40.7 ± 11.8 0 ± 0
8 SSP370 2070_2100 4 ± 5.2 16 ± 5.5 35.5 ± 9.5 13.2 ± 5.2 68.7 ± 13.21 30.9 ± 10.6 0 ± 0
6 SSP585 2030_2060 0 ± 0 13 ± 7 37.6 ± 9 14.1 ± 7.1 64.7 ± 13.43 35 ± 9.4 0 ± 0
9 SSP585 2070_2100 2 ± 2.4 16.6 ± 8 41.2 ± 7.7 13.3 ± 7.7 73.1 ± 13.72 26.6 ± 8.8 0 ± 0
$EUROPE
EUROPE
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0 1.9 ± 1.2 9.6 ± 4.1 4.5 ± 2.8 16 ± 5.11 49 ± 15.3 34.8 ± 20.2
2 historical 1970_2000 0 ± 0 1.5 ± 1 8.9 ± 3.4 3.9 ± 2.1 14.3 ± 4.12 44.9 ± 16.1 40.6 ± 20.4
3 historical 1985_2015 0 ± 0 1.3 ± 1 8.4 ± 2.9 3.8 ± 2 13.5 ± 3.66 46 ± 16.1 40.3 ± 19.9
4 SSP245 2030_2060 0 ± 0 1.9 ± 0.9 10.3 ± 3.6 5.1 ± 1.9 17.3 ± 4.17 52.7 ± 15.4 29.8 ± 19.7
7 SSP245 2070_2100 0 ± 0 2.2 ± 1.1 10.1 ± 4.8 5.1 ± 2.3 17.4 ± 5.44 55.4 ± 14.4 27 ± 18.6
5 SSP370 2030_2060 0 ± 0 2.1 ± 1 9.7 ± 4.7 4.8 ± 2.5 16.6 ± 5.42 52.7 ± 15.7 30.4 ± 20
8 SSP370 2070_2100 0 ± 0 3.1 ± 1.1 11.4 ± 5.9 6.4 ± 3.8 20.9 ± 7.1 52.8 ± 14.3 26 ± 18.6
6 SSP585 2030_2060 0 ± 0 2.3 ± 1.1 10.5 ± 3.7 4.9 ± 2.4 17.7 ± 4.55 55.2 ± 15.8 26.9 ± 19.9
9 SSP585 2070_2100 0 ± 0 3 ± 1.7 11.1 ± 5.4 6.1 ± 2.9 20.2 ± 6.36 57.2 ± 12.8 22.5 ± 17.3
$EUROPE-AFRICA
EUROPE-AFRICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 5.9 ± 4.9 23.6 ± 7 21.3 ± 8.8 8.8 ± 3.3 59.6 ± 12.7 24.9 ± 14.2 14.5 ± 0
2 historical 1970_2000 7 ± 5.5 26.5 ± 7.5 23.3 ± 10.3 9.9 ± 3.2 66.7 ± 14.24 27.5 ± 16.5 4.5 ± 7.5
3 historical 1985_2015 7.1 ± 5.9 26.3 ± 7.1 22.2 ± 7.7 10.3 ± 4.2 65.9 ± 12.73 27.8 ± 14.3 5.2 ± 7
4 SSP245 2030_2060 9.9 ± 6.6 28.9 ± 8.6 26.7 ± 8.9 9.4 ± 3.4 74.9 ± 14.43 23.1 ± 16.2 0.9 ± 0
7 SSP245 2070_2100 11.7 ± 5.7 28.6 ± 8.8 27.2 ± 8.1 9.5 ± 3.1 77 ± 13.61 21.9 ± 13.9 0 ± 0
5 SSP370 2030_2060 10.3 ± 6.3 28.2 ± 8.4 26 ± 8.7 9 ± 2.8 73.5 ± 13.92 24.8 ± 14 0.5 ± 0
8 SSP370 2070_2100 13.2 ± 6.2 27.3 ± 8.6 28.3 ± 9.7 8.7 ± 3.6 77.5 ± 14.81 21.2 ± 16.2 0 ± 0
6 SSP585 2030_2060 11.1 ± 6.4 28.8 ± 8.3 27.5 ± 8.8 8.7 ± 3 76.1 ± 14.01 22.6 ± 14.1 0.2 ± 0
9 SSP585 2070_2100 15.2 ± 5.7 29.1 ± 10.1 27.9 ± 8.8 8 ± 2.8 80.2 ± 14.82 18.6 ± 13.6 0 ± 0
$NORTH-AMERICA
NORTH-AMERICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0 0.4 ± 0.2 5.8 ± 5.6 4.3 ± 2.6 10.5 ± 6.18 35.1 ± 9 54.2 ± 10.9
2 historical 1970_2000 0 ± 0 0.4 ± 0.2 4.3 ± 4.4 4.2 ± 2.4 8.9 ± 5.02 34.5 ± 9.8 56.5 ± 11.7
3 historical 1985_2015 0 ± 0 0.3 ± 0.2 4.6 ± 3.8 4.1 ± 2.2 9 ± 4.4 33.6 ± 10.3 57.2 ± 11.5
4 SSP245 2030_2060 0 ± 0 0.4 ± 0.3 6.5 ± 5.5 4.7 ± 2.6 11.6 ± 6.09 37.5 ± 8.7 50.8 ± 11.1
7 SSP245 2070_2100 0 ± 0 0.5 ± 0.3 6.4 ± 5.1 5.1 ± 2.9 12 ± 5.87 39.7 ± 9 48.1 ± 11
5 SSP370 2030_2060 0 ± 0 0.5 ± 0.4 5.9 ± 4.8 4.6 ± 2.7 11 ± 5.52 37 ± 9.7 51.8 ± 12.2
8 SSP370 2070_2100 0.1 ± 0 0.8 ± 0.7 7.6 ± 4 5.2 ± 2 13.7 ± 4.53 38.4 ± 10.6 47.7 ± 11.2
6 SSP585 2030_2060 0 ± 0 0.5 ± 0.4 6.6 ± 5.2 5.2 ± 2.8 12.3 ± 5.92 38.2 ± 9.1 49.3 ± 11.7
9 SSP585 2070_2100 0.2 ± 0 1 ± 1.4 7 ± 6.1 5.2 ± 2.5 13.4 ± 6.74 42.1 ± 9.3 44.5 ± 11.3
$OCEANIA
OCEANIA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 5.3 ± 0 32 ± 26 36.7 ± 17.1 6.6 ± 3.1 80.6 ± 31.27 18.9 ± 12.5 0.3 ± 0.2
2 historical 1970_2000 5.2 ± 0 31.3 ± 26.4 36.8 ± 17.4 6.8 ± 3.1 80.1 ± 31.77 19.3 ± 12.8 0.5 ± 0.4
3 historical 1985_2015 5.2 ± 0 33.4 ± 26.1 35.5 ± 17.9 6.8 ± 3.2 80.9 ± 31.81 18.4 ± 12.3 0.6 ± 0.3
4 SSP245 2030_2060 10.4 ± 0 35.8 ± 23.5 31.7 ± 15.1 5.7 ± 2.6 83.6 ± 28.05 15.9 ± 11.9 0.3 ± 0.3
7 SSP245 2070_2100 10.5 ± 0 37.9 ± 22.2 29.1 ± 13.1 6.2 ± 3.3 83.7 ± 25.99 16.2 ± 12 0.1 ± 0
5 SSP370 2030_2060 10.7 ± 0 34.8 ± 22.6 33.4 ± 15.3 5.8 ± 2.5 84.7 ± 27.41 14.9 ± 8.1 0.3 ± 0.2
8 SSP370 2070_2100 10.5 ± 10.2 38.1 ± 19.7 28.5 ± 13.8 5.7 ± 3.2 82.8 ± 26.32 16.8 ± 16.3 0.3 ± 0.3
6 SSP585 2030_2060 5.6 ± 3.6 35.6 ± 25.3 33.8 ± 16.7 6.7 ± 3.6 81.7 ± 30.74 17.7 ± 13.8 0.3 ± 0.3
9 SSP585 2070_2100 9.3 ± 6.4 39.4 ± 21.7 29.4 ± 12.9 5.5 ± 2.5 83.6 ± 26.16 16.2 ± 12.5 0 ± 0
$SOUTH-AMERICA
SOUTH-AMERICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0.1 ± 0.1 2.1 ± 1.7 15.2 ± 8.7 7.4 ± 2.6 24.8 ± 9.24 73.2 ± 11.7 1.9 ± 1
2 historical 1970_2000 0.1 ± 0.1 1.8 ± 1.3 13.9 ± 4.3 7.3 ± 2.3 23.1 ± 5.05 74.5 ± 7.3 2.3 ± 2.6
3 historical 1985_2015 0.2 ± 0.1 2 ± 1.2 14.5 ± 5.1 7.5 ± 2.4 24.2 ± 5.76 73.6 ± 7.7 2.1 ± 2.6
4 SSP245 2030_2060 0.2 ± 0.2 2.7 ± 2.2 17.4 ± 9.6 8.8 ± 3.6 29.1 ± 10.49 69.4 ± 13.1 1.3 ± 0.9
7 SSP245 2070_2100 0.2 ± 0.1 3.3 ± 3.5 18.3 ± 9.1 9.3 ± 4.6 31.1 ± 10.78 67.5 ± 14.7 1.3 ± 0.7
5 SSP370 2030_2060 0.2 ± 0.1 2.2 ± 1.7 16.3 ± 5.5 8.8 ± 3.7 27.5 ± 6.84 70.5 ± 8.9 1.9 ± 1.6
8 SSP370 2070_2100 0.2 ± 0.1 3.3 ± 2.5 18.6 ± 7 8.7 ± 3.8 30.8 ± 8.35 67.3 ± 10.8 1.8 ± 1.8
6 SSP585 2030_2060 0.2 ± 0.1 3.1 ± 2.5 18.6 ± 9 9.2 ± 4.4 31.1 ± 10.33 67.6 ± 13.1 1.2 ± 0.8
9 SSP585 2070_2100 0.6 ± 1 4.4 ± 3.8 19.6 ± 9.5 9.7 ± 3.9 34.3 ± 11 64.4 ± 12.6 1.1 ± 0.7

9.3.2 Changes in proportion of aridity categories by continent

kdiff_list
$AFRICA
AFRICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0.3 ± 6.6 0.4 ± 6.5 0.5 ± 5 0.1 ± 2.1 1.3 ± 10.7 0.9 ± 9.4 -2.3 ± 0
2 historical 1970_2000 0 ± 6.6 0 ± 5.9 0 ± 4.8 0 ± 2.1 0 ± 10.4 0 ± 8.8 0 ± 0
3 historical 1985_2015 -0.4 ± 6.7 0.4 ± 6.2 0.3 ± 4.7 -0.2 ± 2.1 0.1 ± 10.5 -0.1 ± 8.8 -0.1 ± 0
4 SSP245 2030_2060 0.1 ± 6.7 1.3 ± 6.5 1.4 ± 4.9 0.2 ± 2.3 3 ± 10.8 -0.8 ± 9.8 -2.3 ± 0
7 SSP245 2070_2100 1 ± 6.6 0.7 ± 6.5 2.1 ± 4.8 0 ± 2.1 3.8 ± 10.7 -1.5 ± 10 -2.3 ± 0
5 SSP370 2030_2060 -0.9 ± 6.9 -0.1 ± 5.5 2.1 ± 3.8 0.6 ± 1.8 1.7 ± 9.8 -0.5 ± 9.1 -1.2 ± 0
8 SSP370 2070_2100 1.1 ± 5.9 0.9 ± 6.4 1.5 ± 4.6 0 ± 2.2 3.5 ± 10.1 -2.4 ± 8.9 -1.1 ± 0
6 SSP585 2030_2060 0.2 ± 6.6 0.9 ± 6.4 1.7 ± 4.9 0.1 ± 2.1 2.9 ± 10.6 -0.7 ± 9.7 -2.3 ± 0
9 SSP585 2070_2100 2.8 ± 11.5 0.6 ± 6.8 1.1 ± 5.4 -0.2 ± 2.1 4.3 ± 14.6 -1.9 ± 10.6 -2.3 ± 0
$ASIA
ASIA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0.1 ± 3.1 0.2 ± 2.6 0.3 ± 2.1 0.2 ± 1.8 0.8 ± 4.9 1.7 ± 11.5 -2.5 ± 13.3
2 historical 1970_2000 0 ± 3 0 ± 2.5 0 ± 2.1 0 ± 1.6 0 ± 4.7 0 ± 12 0 ± 13.6
3 historical 1985_2015 -0.1 ± 3 -0.2 ± 2.4 -0.3 ± 2.1 -0.1 ± 1.6 -0.7 ± 4.7 -0.6 ± 11.4 1.3 ± 13
4 SSP245 2030_2060 -0.3 ± 2.8 0.5 ± 2.4 0.5 ± 2.2 0.5 ± 2 1.2 ± 4.7 3.6 ± 11.7 -4.7 ± 13.3
7 SSP245 2070_2100 -0.1 ± 2.9 0.5 ± 2.5 0.6 ± 2.3 0.4 ± 2 1.4 ± 4.9 5.5 ± 12.1 -6.8 ± 13.4
5 SSP370 2030_2060 -0.4 ± 3.1 0.4 ± 2.3 0 ± 2.3 -0.1 ± 1.7 -0.1 ± 4.8 3.3 ± 11.1 -3.2 ± 13.5
8 SSP370 2070_2100 0.8 ± 2.8 0.5 ± 2.3 0.7 ± 2.1 0 ± 1.6 2 ± 4.5 4.5 ± 12.5 -6.4 ± 13.6
6 SSP585 2030_2060 0 ± 2.9 0.5 ± 2.4 0.6 ± 2.3 0.4 ± 1.9 1.5 ± 4.8 4.7 ± 12 -6.1 ± 13.4
9 SSP585 2070_2100 1 ± 3.8 1.1 ± 3.7 1.5 ± 3.7 0.4 ± 1.9 4 ± 6.7 5.8 ± 14.5 -9.8 ± 13.8
$CENTRAL-AMERICA
CENTRAL-AMERICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0 0 ± 5.5 0.4 ± 9.2 1.4 ± 11.7 1.8 ± 15.8 0.1 ± 13.5 -2 ± 0
2 historical 1970_2000 0 ± 0 0 ± 5.8 0 ± 10.7 0 ± 11 0 ± 16.5 0 ± 17 0 ± 0
3 historical 1985_2015 0 ± 0 0.6 ± 5.9 1 ± 10.3 0.1 ± 10.8 1.7 ± 16 -1.1 ± 16.6 -0.6 ± 0
4 SSP245 2030_2060 0 ± 0 3.5 ± 7.4 6.6 ± 11 -0.8 ± 10 9.3 ± 16.6 -7.3 ± 14.8 -2 ± 0
7 SSP245 2070_2100 0 ± 0 3.9 ± 6.9 8.4 ± 11.4 -0.3 ± 11.1 12 ± 17.3 -10.1 ± 14.8 -2 ± 0
5 SSP370 2030_2060 0 ± 0 4.2 ± 7.5 1.7 ± 10.4 0.3 ± 11.1 6.2 ± 17 -4.2 ± 16.8 -2 ± 0
8 SSP370 2070_2100 4 ± 5.2 7.6 ± 6.9 5.8 ± 12.2 -1.5 ± 9.4 15.9 ± 17.6 -14 ± 16 -2 ± 0
6 SSP585 2030_2060 0 ± 0 4.6 ± 8.1 7.9 ± 11.8 -0.6 ± 10.5 11.9 ± 17.8 -9.9 ± 15.2 -2 ± 0
9 SSP585 2070_2100 2 ± 2.4 8.2 ± 9 11.5 ± 10.8 -1.4 ± 11 20.3 ± 18 -18.3 ± 14.9 -2 ± 0
$EUROPE
EUROPE
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0 0.4 ± 1.6 0.7 ± 5.3 0.6 ± 3.5 1.7 ± 6.6 4.1 ± 22.2 -5.8 ± 28.7
2 historical 1970_2000 0 ± 0 0 ± 1.4 0 ± 4.8 0 ± 3 0 ± 5.8 0 ± 22.8 0 ± 28.8
3 historical 1985_2015 0 ± 0 -0.2 ± 1.4 -0.5 ± 4.5 -0.1 ± 2.9 -0.8 ± 5.5 1.1 ± 22.8 -0.3 ± 28.5
4 SSP245 2030_2060 0 ± 0 0.4 ± 1.3 1.4 ± 5 1.2 ± 2.8 3 ± 5.9 7.8 ± 22.3 -10.8 ± 28.4
7 SSP245 2070_2100 0 ± 0 0.7 ± 1.5 1.2 ± 5.9 1.2 ± 3.1 3.1 ± 6.8 10.5 ± 21.6 -13.6 ± 27.6
5 SSP370 2030_2060 0 ± 0 0.6 ± 1.4 0.8 ± 5.8 0.9 ± 3.3 2.3 ± 6.8 7.8 ± 22.5 -10.2 ± 28.6
8 SSP370 2070_2100 0 ± 0 1.6 ± 1.5 2.5 ± 6.8 2.5 ± 4.3 6.6 ± 8.2 7.9 ± 21.5 -14.6 ± 27.6
6 SSP585 2030_2060 0 ± 0 0.8 ± 1.5 1.6 ± 5 1 ± 3.2 3.4 ± 6.1 10.3 ± 22.6 -13.7 ± 28.5
9 SSP585 2070_2100 0 ± 0 1.5 ± 2 2.2 ± 6.4 2.2 ± 3.6 5.9 ± 7.6 12.3 ± 20.6 -18.1 ± 26.7
$EUROPE-AFRICA
EUROPE-AFRICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 -1.1 ± 7.4 -2.9 ± 10.3 -2 ± 13.5 -1.1 ± 4.6 -7.1 ± 19.1 -2.6 ± 21.8 10 ± 7.5
2 historical 1970_2000 0 ± 7.8 0 ± 10.6 0 ± 14.6 0 ± 4.5 0 ± 20.1 0 ± 23.3 0 ± 10.6
3 historical 1985_2015 0.1 ± 8.1 -0.2 ± 10.3 -1.1 ± 12.9 0.4 ± 5.3 -0.8 ± 19.1 0.3 ± 21.8 0.7 ± 10.3
4 SSP245 2030_2060 2.9 ± 8.6 2.4 ± 11.4 3.4 ± 13.6 -0.5 ± 4.7 8.2 ± 20.3 -4.4 ± 23.1 -3.6 ± 7.5
7 SSP245 2070_2100 4.7 ± 7.9 2.1 ± 11.6 3.9 ± 13.1 -0.4 ± 4.5 10.3 ± 19.7 -5.6 ± 21.6 -4.5 ± 7.5
5 SSP370 2030_2060 3.3 ± 8.4 1.7 ± 11.3 2.7 ± 13.5 -0.9 ± 4.3 6.8 ± 19.9 -2.7 ± 21.6 -4 ± 7.5
8 SSP370 2070_2100 6.2 ± 8.3 0.8 ± 11.4 5 ± 14.1 -1.2 ± 4.8 10.8 ± 20.5 -6.3 ± 23.1 -4.5 ± 7.5
6 SSP585 2030_2060 4.1 ± 8.4 2.3 ± 11.2 4.2 ± 13.5 -1.2 ± 4.4 9.4 ± 20 -4.9 ± 21.7 -4.3 ± 7.5
9 SSP585 2070_2100 8.2 ± 7.9 2.6 ± 12.6 4.6 ± 13.5 -1.9 ± 4.3 13.5 ± 20.6 -8.9 ± 21.4 -4.5 ± 7.5
$NORTH-AMERICA
NORTH-AMERICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0 0 ± 0.3 1.5 ± 7.1 0.1 ± 3.5 1.6 ± 8 0.6 ± 13.3 -2.3 ± 16
2 historical 1970_2000 0 ± 0 0 ± 0.3 0 ± 6.2 0 ± 3.4 0 ± 7.1 0 ± 13.9 0 ± 16.5
3 historical 1985_2015 0 ± 0 -0.1 ± 0.3 0.3 ± 5.8 -0.1 ± 3.3 0.1 ± 6.7 -0.9 ± 14.2 0.7 ± 16.4
4 SSP245 2030_2060 0 ± 0 0 ± 0.4 2.2 ± 7 0.5 ± 3.5 2.7 ± 7.9 3 ± 13.1 -5.7 ± 16.1
7 SSP245 2070_2100 0 ± 0 0.1 ± 0.4 2.1 ± 6.7 0.9 ± 3.8 3.1 ± 7.7 5.2 ± 13.3 -8.4 ± 16.1
5 SSP370 2030_2060 0 ± 0 0.1 ± 0.4 1.6 ± 6.5 0.4 ± 3.6 2.1 ± 7.5 2.5 ± 13.8 -4.7 ± 16.9
8 SSP370 2070_2100 0.1 ± 0 0.4 ± 0.7 3.3 ± 5.9 1 ± 3.1 4.8 ± 6.8 3.9 ± 14.4 -8.8 ± 16.2
6 SSP585 2030_2060 0 ± 0 0.1 ± 0.4 2.3 ± 6.8 1 ± 3.7 3.4 ± 7.8 3.7 ± 13.4 -7.2 ± 16.5
9 SSP585 2070_2100 0.2 ± 0 0.6 ± 1.4 2.7 ± 7.5 1 ± 3.5 4.5 ± 8.4 7.6 ± 13.5 -12 ± 16.3
$OCEANIA
OCEANIA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0.1 ± 0 0.7 ± 37.1 -0.1 ± 24.4 -0.2 ± 4.4 0.5 ± 44.6 -0.4 ± 17.9 -0.2 ± 0.4
2 historical 1970_2000 0 ± 0 0 ± 37.3 0 ± 24.6 0 ± 4.4 0 ± 44.9 0 ± 18.1 0 ± 0.6
3 historical 1985_2015 0 ± 0 2.1 ± 37.1 -1.3 ± 25 0 ± 4.5 0.8 ± 45 -0.9 ± 17.8 0.1 ± 0.5
4 SSP245 2030_2060 5.2 ± 0 4.5 ± 35.3 -5.1 ± 23 -1.1 ± 4 3.5 ± 42.4 -3.4 ± 17.5 -0.2 ± 0.5
7 SSP245 2070_2100 5.3 ± 0 6.6 ± 34.5 -7.7 ± 21.8 -0.6 ± 4.5 3.6 ± 41 -3.1 ± 17.5 -0.4 ± 0.4
5 SSP370 2030_2060 5.5 ± 0 3.5 ± 34.8 -3.4 ± 23.2 -1 ± 4 4.6 ± 42 -4.4 ± 15.1 -0.2 ± 0.4
8 SSP370 2070_2100 5.3 ± 10.2 6.8 ± 32.9 -8.3 ± 22.2 -1.1 ± 4.5 2.7 ± 41.3 -2.5 ± 20.7 -0.2 ± 0.5
6 SSP585 2030_2060 0.4 ± 3.6 4.3 ± 36.6 -3 ± 24.1 -0.1 ± 4.8 1.6 ± 44.2 -1.6 ± 18.8 -0.2 ± 0.5
9 SSP585 2070_2100 4.1 ± 6.4 8.1 ± 34.2 -7.4 ± 21.7 -1.3 ± 4 3.5 ± 41.2 -3.1 ± 17.9 -0.5 ± 0.4
$SOUTH-AMERICA
SOUTH-AMERICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0.1 0.3 ± 2.1 1.3 ± 9.7 0.1 ± 3.5 1.7 ± 10.5 -1.3 ± 13.8 -0.4 ± 2.8
2 historical 1970_2000 0 ± 0.1 0 ± 1.8 0 ± 6.1 0 ± 3.3 0 ± 7.1 0 ± 10.3 0 ± 3.7
3 historical 1985_2015 0.1 ± 0.1 0.2 ± 1.8 0.6 ± 6.7 0.2 ± 3.3 1.1 ± 7.7 -0.9 ± 10.6 -0.2 ± 3.7
4 SSP245 2030_2060 0.1 ± 0.2 0.9 ± 2.6 3.5 ± 10.5 1.5 ± 4.3 6 ± 11.6 -5.1 ± 15 -1 ± 2.8
7 SSP245 2070_2100 0.1 ± 0.1 1.5 ± 3.7 4.4 ± 10.1 2 ± 5.1 8 ± 11.9 -7 ± 16.4 -1 ± 2.7
5 SSP370 2030_2060 0.1 ± 0.1 0.4 ± 2.1 2.4 ± 7 1.5 ± 4.4 4.4 ± 8.5 -4 ± 11.5 -0.4 ± 3.1
8 SSP370 2070_2100 0.1 ± 0.1 1.5 ± 2.8 4.7 ± 8.2 1.4 ± 4.4 7.7 ± 9.8 -7.2 ± 13 -0.5 ± 3.2
6 SSP585 2030_2060 0.1 ± 0.1 1.3 ± 2.8 4.7 ± 10 1.9 ± 5 8 ± 11.5 -6.9 ± 15 -1.1 ± 2.7
9 SSP585 2070_2100 0.5 ± 1 2.6 ± 4 5.7 ± 10.4 2.4 ± 4.5 11.2 ± 12.1 -10.1 ± 14.6 -1.2 ± 2.7

9.3.3 Changes in proportion of aridity categories, by continent, calculated by model


cmip6$cat.AI <- factor(cmip6$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

tab.percent.cont <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, Continent, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() %>% 
  group_by(period, model, Continent, cat.AI) %>% 
  summarise(mmmean = mean(count, na.rm = T), mmsd = sd(count, na.rm = T)) %>%
  mutate(percent = round(mmmean/sum(mmmean)*100, 1), sd.percent = round(mmsd/sum(mmmean)*100,1)) 
write.table(tab.percent.cont, "tab.percent.continent.txt")

tab.cont <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, Continent, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() 

list_ref <- list()

for(i in unique(tab.cont$source)){
  tab.cont.i <- filter(tab.cont, source == i)
  tab.cont.i.ref <- filter(tab.cont.i, period == "1970_2000")
  
  for(j in unique(tab.cont.i$Continent)){
    tab.cont.i.j <- filter(tab.cont.i, Continent == j)
    tab.cont.i.j.ref <- filter(tab.cont.i.ref, Continent == j)
    
    for(k in unique(tab.cont.i$cat.AI)){
      
      index <- paste(i, j, k, sep = "_")
      tab.cont.i.j.k <- filter(tab.cont.i.j, cat.AI == k)
      
      lines <- ifelse(dim(tab.cont.i.j.k)[1] > 0, dim(tab.cont.i.j.k)[1] - 1, 1)
      
      tab.cont.i.j.k.ref <- subset(tab.cont.i.j.ref, cat.AI == k)
      
      if(dim(tab.cont.i.j.k.ref)[1] == 0){
        names1 <- names(tab.cont.i.j.k.ref)
        line1 <- c("1970_2000","historical", i, j, k, 0)
        tab.cont.i.j.k.ref <- rbind(tab.cont.i.j.k.ref, line1) %>%
          setNames(names1)
        
      }else{
        tab.cont.i.j.k.ref <- tab.cont.i.j.k.ref
      }
      
      tab.cont.i.j.k.ref <- tab.cont.i.j.k.ref %>% rbind(tab.cont.i.j.k.ref[rep(1,lines),])
      tab.cont.i.j.k$diff <- tab.cont.i.j.k$count - as.numeric(tab.cont.i.j.k.ref$count)  
      
      list_ref[[index]] <- tab.cont.i.j.k
    }
  }
}


tab.cont.diff <- bind_rows(list_ref, .id = "column_label") %>%
  group_by(period, model, Continent, cat.AI) %>% 
  summarise(diff.mean = mean(diff, na.rm = T), diff.sd = sd(diff, na.rm = T))
write.table(tab.cont.diff, "tab.cont.diff.txt")
tab.cont.diff <- read.table("tab.cont.diff.txt")

tab_list_percent <- list()
tab_list_sd <- list()

for(i in unique(tab.cont.diff$Continent)){
  tab.i <- subset(tab.cont.diff, Continent == i)
  
  df.percent <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "diff.mean")
  
  df.sd <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "diff.sd")
  # sd for a sum: square root of sum of squared sd
  
  drycats <- which(names(df.percent) %in% c("Hyperarid","Arid","Semi-arid","Dry subhumid"))
  
  df.percent <- df.percent %>% mutate("Sum drylands" = rowSums(.[drycats], na.rm = T))
  df.sd <- df.sd %>% mutate("Sum drylands" = sqrt(rowSums((.[drycats])^2, na.rm = T)),2)
  
  missing <- setdiff(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"), names(df.percent))
  
  df.percent[, missing] <- 0
  df.sd[,missing] <- 0
  
  df.percent <- df.percent %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.percent[is.na(df.percent)] <- 0
  df.percent[,c(3:9)] <- round(df.percent[,c(3:9)], 1)
  
  tab_list_percent[[i]] <- df.percent
  
  df.sd <- df.sd %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.sd[is.na(df.sd)] <- 0
  df.sd[ , c(3:9)] <- round(df.sd[, c(3:9)], 1)
  
  tab_list_sd[[i]] <- df.sd
}

k_list <- list()

for(i in names(tab_list_percent)){
  
  df <- tab_list_percent[[i]] %>% select(c("model", "period"))
  
  df$Hyperarid <- paste(tab_list_percent[[i]]$Hyperarid, tab_list_sd[[i]]$Hyperarid, sep = " ± ")
  df$Arid <- paste(tab_list_percent[[i]]$Arid, tab_list_sd[[i]]$Arid, sep = " ± ")
  df$'Semi-Arid' <- paste(tab_list_percent[[i]][,5], tab_list_sd[[i]][,5], sep = " ± ")
  df$'Dry subhumid' <- paste(tab_list_percent[[i]][,6], tab_list_sd[[i]][,6], sep = " ± ")
  df$'Sum drylands' <- paste(tab_list_percent[[i]][,7], tab_list_sd[[i]][,7], sep = " ± ")
  df$Humid <- paste(tab_list_percent[[i]]$Humid, tab_list_sd[[i]]$Humid, sep = " ± ")
  df$Cold <- paste(tab_list_percent[[i]]$Cold, tab_list_sd[[i]]$Cold, sep = " ± ")
  
  k_list[[i]] <- kable(df[order(df$model, decreasing = F)], caption = i, digits = 2) %>% kable_styling(bootstrap_options = "bordered") %>%
    column_spec(8, italic = T, background = "#DAE2F7") %>% row_spec(2, bold = T) 
  
}

k_list
$AFRICA
AFRICA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 -5.8 ± 24.1 3 ± 22.7 1.5 ± 48.9 4.8 ± 34 2.8 ± 38.7 -0.4 ± 12 0 ± 0
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 -10.5 ± 19.5 11 ± 21.6 2.8 ± 33.1 5.9 ± 11.9 -2.7 ± 15.5 -3.5 ± 10.1 -2 ± 0
SSP245 2030_2060 -10.5 ± 41 24.3 ± 29.7 42.2 ± 66.6 25.8 ± 40.3 -37.8 ± 53.2 2.5 ± 15.5 0 ± 0
SSP370 2030_2060 -29.5 ± 88.4 -5.2 ± 60.9 23.7 ± 145 45.5 ± 86.1 -21.4 ± 43.4 12.9 ± 45.6 -30 ± 0
SSP585 2030_2060 -7.8 ± 50.9 13.6 ± 39.6 39 ± 82.2 33.2 ± 47.6 -34.6 ± 52.7 0 ± 18.4 0 ± 0
SSP245 2070_2100 10.1 ± 66.8 9 ± 46.2 58.8 ± 99 42.1 ± 51.9 -54.5 ± 61.7 -2.3 ± 22.7 0 ± 0
SSP370 2070_2100 19.5 ± 66.3 18.5 ± 58.2 69.2 ± 103.5 31.1 ± 46.6 -67.1 ± 52.3 0.2 ± 27.4 -28 ± 0
SSP585 2070_2100 45.6 ± 209.6 9.2 ± 79.8 67.5 ± 240.3 16.5 ± 83.8 -62.8 ± 98 -3.8 ± 20.8 0 ± 0
$ASIA
ASIA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 4.3 ± 28 8.8 ± 41.5 40.2 ± 88.6 11.3 ± 61.5 90.2 ± 176.8 15.8 ± 39.6 -130.5 ± 287.7
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 -5.2 ± 20 -12.6 ± 42 -36.4 ± 57.8 -17 ± 29.8 -33.6 ± 169.1 -1.6 ± 17.2 70 ± 234.1
SSP245 2030_2060 -15.7 ± 38.1 24.2 ± 40.5 57.8 ± 123.7 22.2 ± 88.9 192.7 ± 202.9 27 ± 65.5 -250.5 ± 375.8
SSP370 2030_2060 -23.8 ± 61.4 21.6 ± 45.7 -4.5 ± 98.3 0 ± 55.9 174.2 ± 120.1 -2.3 ± 26 -169.7 ± 100.4
SSP585 2030_2060 -1.8 ± 30.8 23.8 ± 34.1 72.9 ± 115.9 27.1 ± 91.2 250 ± 201.3 23.8 ± 54.9 -322.9 ± 358.7
SSP245 2070_2100 -8.5 ± 31.9 23.6 ± 37.9 68.1 ± 133.8 30.1 ± 101.2 294.5 ± 208.7 22.8 ± 72.2 -362.5 ± 381.8
SSP370 2070_2100 40.1 ± 64.3 25.3 ± 63.8 101.2 ± 108.9 32.7 ± 53 241.2 ± 179.9 3.1 ± 29 -342.3 ± 156
SSP585 2070_2100 55.1 ± 178.4 59.5 ± 146 215 ± 313.2 77.9 ± 202.8 324.1 ± 262.6 22.5 ± 61.9 -539.1 ± 416.4
$CENTRAL-AMERICA
CENTRAL-AMERICA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 0 ± 0 -0.3 ± 4.7 2.5 ± 21.5 -0.5 ± 18 -2 ± 28.2 3.3 ± 10.8 0 ± 0
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 0 ± 0 1.7 ± 4 4.2 ± 8.6 2.5 ± 6.3 -4 ± 5 -0.1 ± 4.3 -2 ± 0
SSP245 2030_2060 0 ± 0 9.8 ± 11.1 24 ± 29.7 17.3 ± 25.7 -23.5 ± 26.8 -3.2 ± 10 0 ± 0
SSP370 2030_2060 0 ± 0 11.8 ± 12.4 15.1 ± 31.2 3.3 ± 24 -14.6 ± 37.3 0 ± 15.6 0 ± 0
SSP585 2030_2060 0 ± 0 12.8 ± 14 31.5 ± 36.2 21.1 ± 29.1 -31 ± 25.5 -2.5 ± 16.4 0 ± 0
SSP245 2070_2100 0 ± 0 11 ± 10.3 31.9 ± 36.4 22.6 ± 29.2 -31.5 ± 27.3 -1.7 ± 19 0 ± 0
SSP370 2070_2100 12 ± 15.6 23.2 ± 15.1 50.2 ± 39 18.7 ± 29.5 -39.6 ± 23.8 -3.7 ± 13.4 0 ± 0
SSP585 2070_2100 6 ± 7.1 25.2 ± 24.9 60.6 ± 42.8 35.6 ± 28.6 -55.1 ± 42.4 -6.2 ± 18.6 0 ± 0
$EUROPE
EUROPE
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 0 ± 0 4.7 ± 11.2 20.7 ± 20.7 8.2 ± 13.3 50.3 ± 106.6 7.8 ± 11.3 -70 ± 115
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 0 ± 0 -2.2 ± 9.8 -9.2 ± 27.8 -5.8 ± 24.6 12.9 ± 48.5 -1.2 ± 8.5 -4.3 ± 73.1
SSP245 2030_2060 0 ± 0 4.7 ± 7.4 33.7 ± 18.9 17.5 ± 14.6 97.5 ± 126 11.5 ± 9.6 -128.5 ± 146.3
SSP370 2030_2060 0 ± 0 7.2 ± 8.5 35.2 ± 16.1 17.6 ± 11.2 90.5 ± 75.1 10.4 ± 7.8 -125.2 ± 73.5
SSP585 2030_2060 0 ± 0 9.8 ± 8.2 41.1 ± 18.4 19.1 ± 14.2 125.9 ± 134.4 12.2 ± 8.4 -164.8 ± 144.4
SSP245 2070_2100 0 ± 0 8.6 ± 11.2 44.2 ± 24.7 21.6 ± 19.9 123.1 ± 120.8 14 ± 9.3 -166.6 ± 139.5
SSP370 2070_2100 0 ± 0 19.2 ± 13.7 86.9 ± 44.6 37.7 ± 28 92.2 ± 110 30 ± 31.8 -177.7 ± 82.2
SSP585 2070_2100 0 ± 0 18.3 ± 19 81.8 ± 39.7 37.4 ± 24.7 150.7 ± 131 26.1 ± 24.6 -230.9 ± 148.4
$EUROPE-AFRICA
EUROPE-AFRICA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 -2.1 ± 5.2 -1.2 ± 9 -1.8 ± 19.9 1.8 ± 14.9 1 ± 20.7 -0.4 ± 8.2 13 ± 0
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 0.5 ± 5.8 0.2 ± 8.2 -1.4 ± 24.4 -4.2 ± 17.9 2.5 ± 16.7 2 ± 13.2 -6.5 ± 9.2
SSP245 2030_2060 11.2 ± 8.1 7 ± 10 26.5 ± 25.8 11.6 ± 19.8 -22.2 ± 35.5 -3.3 ± 10.1 -54 ± 0
SSP370 2030_2060 14 ± 6.8 6.1 ± 9.9 26.1 ± 20.3 10.5 ± 14.3 -23.1 ± 20.7 -4.5 ± 8 -56 ± 0
SSP585 2030_2060 16.9 ± 9.9 7.8 ± 10.2 35.5 ± 25.8 16.5 ± 17.2 -33.2 ± 32 -5.8 ± 13 -57 ± 0
SSP245 2070_2100 19.4 ± 8.2 6.8 ± 11.4 38.6 ± 29 14.9 ± 23.1 -36.5 ± 35.9 -2.5 ± 10.7 0 ± 0
SSP370 2070_2100 25.3 ± 16.7 -0.9 ± 23.1 35.5 ± 41.9 17.7 ± 22.1 -30.9 ± 28.5 -6.5 ± 21.3 0 ± 0
SSP585 2070_2100 32.4 ± 12.3 11.3 ± 19.9 56.3 ± 36.9 21.2 ± 22.1 -55.5 ± 27.4 -8.7 ± 18.1 0 ± 0
$NORTH-AMERICA
NORTH-AMERICA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 0 ± 0 0.4 ± 3.1 37 ± 71.9 31.8 ± 66.3 21.7 ± 91.3 4.9 ± 27.7 -56.1 ± 174
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 0 ± 0 0.3 ± 4.2 -5.4 ± 48.8 -3 ± 32.2 -19.8 ± 70 -2.7 ± 36.4 25.1 ± 114
SSP245 2030_2060 0 ± 0 3.4 ± 2.9 77.5 ± 109.8 60 ± 104.2 82.2 ± 94.5 14.2 ± 34.7 -159.2 ± 214.3
SSP370 2030_2060 0 ± 0 5.5 ± 6.9 62 ± 56.5 44.5 ± 45.2 70.5 ± 48.5 12.1 ± 33.2 -131.6 ± 77
SSP585 2030_2060 0 ± 0 7.1 ± 7.9 99 ± 102.3 62.1 ± 91 103.4 ± 97.9 29.8 ± 46.2 -201.3 ± 203.2
SSP245 2070_2100 0 ± 0 4.7 ± 3.1 87.1 ± 99.3 56.9 ± 86 145.9 ± 122.4 25.5 ± 49.7 -231.9 ± 215.6
SSP370 2070_2100 4 ± 0 15.7 ± 20 139.6 ± 80.2 91.5 ± 59.6 108 ± 110.3 28.4 ± 49.8 -243.9 ± 106
SSP585 2070_2100 5 ± 0 20.2 ± 38.4 147.1 ± 145.7 91.2 ± 136 198.8 ± 99.8 30.8 ± 35.4 -339.6 ± 220.2
$OCEANIA
OCEANIA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 1 ± 0 6 ± 42.8 4.5 ± 57.7 -1.2 ± 38.1 -3 ± 12.1 -1.2 ± 6.8 0.2 ± 1.2
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 1 ± 0 5.1 ± 32.9 5.6 ± 40.3 -2.2 ± 21.3 -3 ± 7.9 1.8 ± 9.4 -0.2 ± 1
SSP245 2030_2060 45 ± 0 50.2 ± 62.1 62.8 ± 74.9 -26.5 ± 39.9 -19.8 ± 18.3 -5.8 ± 12.4 -1.5 ± 0.7
SSP370 2030_2060 47 ± 0 42.4 ± 50.3 71.7 ± 61.2 -12.4 ± 33.3 -27 ± 61.2 -5.3 ± 10.3 -3.3 ± 3.2
SSP585 2030_2060 23.5 ± 0.7 33.8 ± 55.9 33.3 ± 71.6 -23.2 ± 43.1 -12 ± 20.1 -0.8 ± 12.2 -1.5 ± 0.7
SSP245 2070_2100 45 ± 0 66.6 ± 81.5 60.7 ± 101.2 -48.5 ± 58.5 -17.4 ± 18.8 -2.5 ± 13 -5 ± 0
SSP370 2070_2100 70.7 ± 84 63.2 ± 110.8 70 ± 153.2 -57 ± 62.9 -14.2 ± 57.5 -6.8 ± 13.5 -5 ± 0
SSP585 2070_2100 60.3 ± 48.7 69.8 ± 126.8 71.2 ± 164.5 -50 ± 91.8 -23.9 ± 22.1 -9 ± 14.2 0 ± 0
$SOUTH-AMERICA
SOUTH-AMERICA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 -0.1 ± 0.4 5.1 ± 13.9 26.3 ± 86.6 20.1 ± 83.4 -20.6 ± 79.5 1.3 ± 18.8 -6.8 ± 35.1
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 0.5 ± 1.2 2.7 ± 5.9 16.7 ± 19.2 10.5 ± 15.7 -14.2 ± 13.2 3 ± 9.3 -2.5 ± 5.8
SSP245 2030_2060 1.6 ± 1.1 13.2 ± 26.6 93.8 ± 100 55.6 ± 93.1 -78.5 ± 97.2 23.5 ± 24.8 -18.8 ± 42
SSP370 2030_2060 1.2 ± 0.9 5.9 ± 19.7 69.4 ± 52.4 38.2 ± 31.9 -59.5 ± 55.6 24.1 ± 36.6 -12.9 ± 18.7
SSP585 2030_2060 1.5 ± 1.1 19.2 ± 29.9 123.7 ± 97.4 72.8 ± 84.5 -107.5 ± 99.4 30.2 ± 38.1 -20.1 ± 41.2
SSP245 2070_2100 1.4 ± 1 22.6 ± 49.9 124 ± 111.5 68.9 ± 90.8 -107 ± 130.1 31 ± 41.3 -22.8 ± 43.7
SSP370 2070_2100 2 ± 2.3 22.9 ± 23.8 121.2 ± 92.7 73.5 ± 82.1 -109.5 ± 109.3 22.7 ± 35.9 -15 ± 15.1
SSP585 2070_2100 7.4 ± 15.3 40.5 ± 51.6 182.8 ± 140.1 96.9 ± 122.8 -160.2 ± 133.1 38 ± 40.7 -29.4 ± 47.8

10 Changes in proportion of aridity categories, figure

10.1 Colplot by continent

tab.flow <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, Continent, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model, Continent) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL)

tab.flow$Continent <- str_replace(tab.flow$Continent, "-", "\n")

df245 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP245")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b245 <- ggplot(data = df245, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  facet_grid(rows = vars(Continent), switch = "y")+
  scale_y_continuous(position = "right")+
  labs(title = "SSP 2-4.5", y = "%", x = "")+
  theme_minimal()

df370 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP370")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b370 <- ggplot(data = df370, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  facet_grid(rows = vars(Continent), switch = "y")+
  scale_y_continuous(position = "right")+  labs(title = "SSP 3-7.0", y = "%", x = "")+
  theme_minimal()


df585 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP585")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b585 <- ggplot(data = df585, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  facet_grid(rows = vars(Continent), switch = "y")+
  scale_y_continuous(position = "right")+  labs(title = "SSP 5-8.5", y = "%", x = "")+
  theme_minimal()


ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)

11 Evolution of aridity index by Continent and period

11.1 Map with labels


tab.percent.region <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period != "1985_2015") %>%
  group_by(period, model, Continent, Name) %>% select(c("diff.AI.mean.percent")) %>%
  summarise(AI.diff.mean = mean(diff.AI.mean.percent, na.rm = T), AI.diff.sd= sd(diff.AI.mean.percent, na.rm = T)) %>%
  ungroup()

ipcc_regions.sf <- st_as_sf(ipcc_regions)

percent.region.sf <- merge(ipcc_regions.sf, tab.percent.region, all.y = T)
pbreaks <- c(-10, 0, 10, 20, 30, 40)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "azure","cornsilk", "#f4a582")
#colscale <- c("#c6dbef","white","#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#993344")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot(data = subset(percent.region.sf, period == j & model == i)) + geom_sf(aes(fill = AI.diff.mean))+
          binned_scale(aesthetics = "fill", breaks = pbreaks, palette = function(x) rev(colscale))+
      borders(colour = "grey60")+
      geom_sf_label(aes(label = round(AI.diff.mean,1)),size = 2) +
      labs(title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

11.2 Point figure

11.2.0.1 By Continent

cmip6s$diff.AI.mean.percent <- with(cmip6s, (AI.mean-AI.ref.mean)/AI.ref.mean*100)

tab.percent.continent <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period != "1985_2015") %>%
  group_by(period, model, Continent) %>% select(c("diff.AI.mean.percent")) %>%
  summarise(AI.diff.mean = mean(diff.AI.mean.percent, na.rm = T), AI.diff.sd= sd(diff.AI.mean.percent, na.rm = T), AI.diff.cv = AI.diff.sd/abs(AI.diff.mean)) %>%
  ungroup()

df <- tab.percent.continent %>% filter(period != "1850_1880" & Continent != "POLAR") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  

dfref$Continent <- str_replace(dfref$Continent, "-", "\n")
dfref$Continent <- factor(dfref$Continent, levels = c("NORTH\nAMERICA", "CENTRAL\nAMERICA", "SOUTH\nAMERICA", "EUROPE", "MED\nBASIN", "AFRICA", "ASIA", "OCEANIA"))


dfref$model[which(dfref$model == "SSP245")] <- "SSP 2-4.5"
dfref$model[which(dfref$model == "SSP370")] <- "SSP 3-7.0"
dfref$model[which(dfref$model == "SSP585")] <- "SSP 5-8.5"

dfref$period[which(dfref$period == "1970_2000")] <- "1970-2000"
dfref$period[which(dfref$period == "2030_2060")] <- "2030-2060"
dfref$period[which(dfref$period == "2070_2100")] <- "2070-2100"
    
ggplot(dfref, aes(x = period, y = AI.diff.mean))+
  geom_point(aes(size = AI.diff.cv, col = AI.diff.cv))+
  geom_line(aes(group = Continent), size = 0.1)+
  facet_grid(rows = vars(model), cols = vars(Continent))+
  scale_size_continuous(breaks = c(0,1,5,10,100), transform = "pseudo_log")+
  scale_color_gradient(low = "gray", high = "#D73027", breaks = c(0,1,5,10,100), transform = "pseudo_log")+
  labs(y = "Difference in aridity index compared to 1970-2000, %", x = "")+
  guides(color = guide_legend("Coefficient of variance"), size = guide_legend("Coefficient of variance"))+
  ylim(-40,25)+
  theme_minimal()+
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 90))

knitr::kable(dfref) %>% kable_styling(bootstrap_options = "bordered")
period model Continent AI.diff.mean AI.diff.sd AI.diff.cv
2030-2060 SSP 2-4.5 AFRICA -1.0179744 12.876170 12.6488143
2030-2060 SSP 2-4.5 ASIA -3.4101720 8.068910 2.3661299
2030-2060 SSP 2-4.5 CENTRAL AMERICA
 -16.643101
BASIN
 -18.634281
AMERICA
  -4.100262
AMERICA
 -10.531839
AMERICA
 -13.424730
BASIN
 -21.449362
AMERICA
  22.197054
AMERICA
 -12.712817
AMERICA
 -20.850032
BASIN
 -25.043787
AMERICA
   0.591772
AMERICA
 -13.182876
AMERICA
 -21.874522
BASIN
 -27.588601
AMERICA
   3.130155
AMERICA
 -12.760628
AMERICA
 -26.142999
BASIN
 -25.843929
AMERICA
 -13.993980
AMERICA
 -13.994989
AMERICA
 -32.095687
BASIN
 -38.314460
AMERICA
  -6.192770
AMERICA
 -19.299671
AMERICA
   0.000000
BASIN
   0.000000
AMERICA
   0.000000
AMERICA
   0.000000
AMERICA
   0.000000
BASIN
   0.000000
AMERICA
   0.000000
AMERICA
   0.000000
AMERICA
   0.000000
BASIN
   0.000000
AMERICA
   0.000000
AMERICA
   0.000000

11.2.0.2 By subregion

names(tab.percent.continent) <- c("period","model","Continent","AI.diff.mean.cont","AI.diff.sd.cont","AI.diff.cv.cont")

tab.percent.subregion <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC","POLAR") & lat > -55 & period != "1985_2015") %>%
  group_by(period, model, Name, Continent) %>% select(c("diff.AI.mean.percent")) %>%
  summarise(AI.diff.mean = mean(diff.AI.mean.percent, na.rm = T), AI.diff.sd= sd(diff.AI.mean.percent, na.rm = T), AI.diff.cv = AI.diff.sd/AI.diff.mean) %>%
  ungroup()

tab.percent.subregion <- merge(tab.percent.subregion, tab.percent.continent, by = c("period","model","Continent"), all.x = T)
  
df <- tab.percent.subregion %>% filter(period != "1850_1880") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  

dfref$Name <- str_replace(dfref$Name, "-", "\n") 

dfref$model[which(dfref$model == "SSP245")] <- "SSP 2-4.5"
dfref$model[which(dfref$model == "SSP370")] <- "SSP 3-7.0"
dfref$model[which(dfref$model == "SSP585")] <- "SSP 5-8.5"

dfref$period[which(dfref$period == "1970_2000")] <- "1970-2000"
dfref$period[which(dfref$period == "2030_2060")] <- "2030-2060"
dfref$period[which(dfref$period == "2070_2100")] <- "2070-2100"

for(i in unique(dfref$Continent)){
df <- filter(dfref, Continent == i)
g <- ggplot(df, aes(x = period))+
  geom_point(aes(y = AI.diff.mean, size = AI.diff.cv, col = AI.diff.cv))+
  geom_line(aes(y = AI.diff.mean, group = Name), size = 0.1)+
  geom_point(aes(y = AI.diff.mean.cont, size = AI.diff.cv.cont), col = "gray", alpha = 0.5)+
  geom_line(aes(y = AI.diff.mean.cont, group = Name), size = 0.1, alpha = 0.5)+
  facet_grid(cols = vars(model), rows = vars(Name), switch = "both")+
  scale_size_continuous(breaks = c(0,1,5,10,100), transform = "pseudo_log")+
 scale_color_gradient(low = "#ABD9E9", high = "#D73027", breaks = c(0,1,5,10,100), transform = "pseudo_log")+
  labs(y = "Difference in aridity index compared to 1970-2000, %", x = "", title = i)+
  guides(color = guide_legend("Standard deviation"), size = guide_legend("Standard deviation"))+
  theme_minimal()+
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 90))+  
  scale_y_continuous(position = "right", limits = c(-40,30))
print(g)
}

11.3 Change AI in boxplot

tab.percent.continent <- cmip6s %>% 
  subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period %in% c("1970_2000","2030_2060","2070_2100"))

df <- tab.percent.continent %>% filter(Continent != "POLAR") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  

dfref$Continent <- str_replace(dfref$Continent, "-", "\n")

ggplot(dfref, aes(x = period, y = diff.AI.mean.percent))+
  geom_boxplot(aes(col = Continent), outliers = F)+
  facet_grid(rows = vars(model), cols = vars(Continent), space = "free_y")+
  scale_color_viridis_d(end = 0.8)+
  labs(y = "Difference in aridity index compared to 1970-2000, %", x = "")+
  guides(color = F)+
  theme_minimal()+
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 90))

11.4 Point figure precipitation and temperature

cmip6s$pr.percent <- with(cmip6s, diff.pr.mean/spr.ref.mean*100)

tab.percent.continent <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period != "1985_2015") %>%
  group_by(period, model, Continent) %>% select(c("diff.t.mean", "pr.percent")) %>%
  summarise(t.diff.mean = mean(diff.t.mean, na.rm = T), t.diff.sd= sd(diff.t.mean, na.rm = T), t.diff.cv = t.diff.sd/abs(t.diff.mean),
            pr.diff.mean = mean(pr.percent, na.rm = T), pr.diff.sd= sd(pr.percent, na.rm = T), pr.diff.cv = pr.diff.sd/abs(pr.diff.mean)) %>%
  ungroup()


df <- tab.percent.continent %>% filter(period != "1850_1880" & Continent != "POLAR") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  

dfref$Continent <- str_replace(dfref$Continent, "-", "\n")
dfref$Continent <- factor(dfref$Continent, levels = c("NORTH\nAMERICA", "CENTRAL\nAMERICA", "SOUTH\nAMERICA", "EUROPE", "MED\nBASIN", "AFRICA", "ASIA", "OCEANIA"))


dfref$model[which(dfref$model == "SSP245")] <- "SSP 2-4.5"
dfref$model[which(dfref$model == "SSP370")] <- "SSP 3-7.0"
dfref$model[which(dfref$model == "SSP585")] <- "SSP 5-8.5"

dfref$period[which(dfref$period == "1970_2000")] <- "1970-2000"
dfref$period[which(dfref$period == "2030_2060")] <- "2030-2060"
dfref$period[which(dfref$period == "2070_2100")] <- "2070-2100"
    
ggplot(dfref, aes(x = period))+
  geom_point(aes(y = t.diff.mean*3, size = t.diff.cv, col = "Temperature, °C"), alpha = 0.5)+
  geom_line(aes(y = t.diff.mean*3, group = Continent, col = "Temperature, °C"), size = 0.1)+
  
  geom_point(aes(y = pr.diff.mean, size = pr.diff.cv, col = "Precipitation, %"), alpha = 0.5)+
  geom_line(aes(y = pr.diff.mean, group = Continent, col = "Precipitation, %"), size = 0.1)+
  scale_size_continuous(breaks = c(0,1,5,10,20), transform = "pseudo_log")+
  scale_color_manual(values = c("#3182bd", "#d6604d"))+
  facet_grid(rows = vars(model), cols = vars(Continent))+
  labs(y = "Precipitation anomalies, %", x = "", size = "Coefficient of variance", color = "")+
  ylim(-45,25)+
  theme_minimal()+
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 90))+
  
  scale_y_continuous(sec.axis = sec_axis(trans = ~./3, name = "Temperature anomalies, °C"))

knitr::kable(dfref) %>% kable_styling(bootstrap_options = "bordered")
period model Continent t.diff.mean t.diff.sd t.diff.cv pr.diff.mean pr.diff.sd pr.diff.cv
2030-2060 SSP 2-4.5 AFRICA 1.631458 0.3166157 0.1940692 7.7904137 12.236848 1.5707572
2030-2060 SSP 2-4.5 ASIA 2.338152 0.6685291 0.2859220 9.6790737 6.609489 0.6828638
2030-2060 SSP 2-4.5 CENTRAL AMERICA
   1.68272
BASIN
   1.98698
AMERICA
   2.62947
AMERICA
   1.45728
AMERICA
   2.03568
BASIN
   2.42553
AMERICA
   3.08982
AMERICA
   1.79399
AMERICA
   2.29305
BASIN
   2.66748
AMERICA
   3.57561
AMERICA
   2.05472
AMERICA
   2.47160
BASIN
   2.89067
AMERICA
   4.00689
AMERICA
   2.26903
AMERICA
   3.47010
BASIN
   4.02250
AMERICA
   5.38537
AMERICA
   3.22043
AMERICA
   4.15576
BASIN
   4.84252
AMERICA
   6.36060
AMERICA
   3.92026
AMERICA
   0.00000
BASIN
   0.00000
AMERICA
   0.00000
AMERICA
   0.00000
AMERICA
   0.00000
BASIN
   0.00000
AMERICA
   0.00000
AMERICA
   0.00000
AMERICA
   0.00000
BASIN
   0.00000
AMERICA
   0.00000
AMERICA
   0.00000

tab.percent.continent <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period != "1985_2015") %>%
  group_by(period, model, Name) %>% select(c("diff.t.mean", "pr.percent")) %>%
  summarise(t.diff.mean = mean(diff.t.mean, na.rm = T), t.diff.sd= sd(diff.t.mean, na.rm = T),
            pr.diff.mean = mean(pr.percent, na.rm = T), pr.diff.sd= sd(pr.percent, na.rm = T)) %>%
  ungroup()


df <- tab.percent.continent %>% filter(period != "1850_1880") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  

dfref$Name <- str_replace(dfref$Name, "-", "\n") 

ggplot(dfref, aes(x = period))+
  geom_point(aes(y = t.diff.mean*3, size = t.diff.sd, col = "Temperature, °C"), alpha = 0.5)+
  geom_line(aes(y = t.diff.mean*3, group = Name, col = "Temperature, °C"), size = 0.1)+
  
  geom_point(aes(y = pr.diff.mean, size = pr.diff.sd, col = "Precipitation, %"), alpha = 0.5)+
  geom_line(aes(y = pr.diff.mean, group = Name, col = "Precipitation, %"), size = 0.1)+
  
  scale_color_manual(values = c("#3182bd", "#d6604d"))+
  facet_grid(cols = vars(model), rows = vars(Name), switch = "both")+
  labs(y = "Precipitation anomalies, %", x = "", size = "Standard deviation", color = "")+
  ylim(-45,25)+
  theme_minimal()+
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 90))+
  
  scale_y_continuous(sec.axis = sec_axis(trans = ~./3, name = "Temperature anomalies, °C"))

12 Changes of precipitation and temperature over time and scenarios

12.1 Temperature

12.1.1 Absolute temperature

quantile(cmip6s$t.mean, probs = seq(0,1,0.1)) 
##          0%         10%         20%         30%         40%         50% 
## -53.4587313 -42.5828261 -28.4644864 -17.7502317  -7.2080554  -0.2580962 
##         60%         70%         80%         90%        100% 
##   7.1198406  15.3318920  22.8072968  26.4971129  34.3727997

t.breaks <- c(-50,-40,-30,-20,-10,0,10,20,30)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")
map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = t.mean))+
    borders(colour = "grey60")+
    binned_scale(aesthetics = "fill", breaks = t.breaks, palette = function(x) colscale,
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "Mean annual temperature, °C")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = t.mean))+
      borders(colour = "grey60")+
      binned_scale(aesthetics = "fill", breaks = t.breaks, palette = function(x) colscale,
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Mean annual temperature, °C", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}



ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")


#labels = c("2030-2060", "2070-2100"

12.1.2 Temperature anomalies

quantile(cmip6s$diff.t.mean, probs = seq(0,1,0.1)) 
##         0%        10%        20%        30%        40%        50%        60% 
## -3.5775507 -0.5616741  0.0000000  1.1670755  1.8544034  2.2566573  2.6600679 
##        70%        80%        90%       100% 
##  3.1444949  3.7454555  4.5875606 10.7627615

t.breaks <- c(-3,-0.5,0,1.5,2,3,4,5,10)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")
map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = diff.t.mean))+
    borders(colour = "grey60")+
   scale_fill_gradient2(low = "#08519c",mid = "white", high = "#b2182b", midpoint = 0, transform = "pseudo_log", breaks = c(-2, -1, 0, 1, 2, 3, 4))+
      labs(fill = "Precipitation anomalies in %", title = i)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = diff.t.mean))+
      borders(colour = "grey60")+
   scale_fill_gradient2(low = "#08519c",mid = "white", high = "#b2182b", midpoint = 0, transform = "pseudo_log", breaks = c(-2, -1, 0, 1, 2, 3, 4))+
      labs(fill = "Temperature anomalies in °C", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

### Binned anomalies

YOR <- brewer.pal(n = 11, name ="YlOrRd")
coltemp <- c("[0,1]" = YOR[1], "(1,2]" = YOR[2], "(2,3]" = YOR[3], "(3,4]" = YOR[4], "(4,5]" = YOR[5], "(5,6]" = YOR[6], "(6,7]" = YOR[7], "(7,8]" = YOR[8], "(8,9]" = YOR[9], "(9,10]" = "grey20")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot(data = subset(cmip6s, period == j & model == i), 
                aes(x=lon, y = lat,  fill = cut(diff.t.mean, breaks = c(0,1,2,3,4,5,6,7,8,9,10), include.lowest = T))) + 
      geom_tile(show.legend = T)+
      borders(colour = "grey60")+
      scale_fill_manual(values = coltemp, drop = F)+
      labs(fill = "Temperature\nanomalies in °C", title = index)+
      theme_void()+ylim(-55,90)+
     theme(legend.position = "bottom", legend.text =  element_text(angle = 0))+
    guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

12.2 Precipitation

12.2.1 Absolute values of precipitations

quantile(cmip6s$pr.mean, probs = seq(0,1,0.1)) 
##          0%         10%         20%         30%         40%         50% 
##    8.579696   55.000319  104.132310  200.852004  337.042710  463.946055 
##         60%         70%         80%         90%        100% 
##  579.741473  712.099280  959.490392 1431.315602 6863.723388

p.breaks <- c(10,50,100,200,300,400,500,600,700)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")
map_list <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = pr.mean))+
    borders(colour = "grey60")+
    binned_scale(aesthetics = "fill", breaks = p.breaks, palette = function(x) rev(colscale),
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "Annual mean\nprecipitation, mm")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = pr.mean))+
      borders(colour = "grey60")+
      binned_scale(aesthetics = "fill", breaks = p.breaks, palette = function(x) rev(colscale),
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Annual mean\nprecipitation, mm", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

12.2.2 Precipitation anomalies

12.2.2.1 Gradient anomalies

cmip6s$pr.percent <- with(cmip6s, diff.pr.mean/spr.ref.mean*100)

quantile(cmip6s$pr.percent, probs = seq(0,1,0.1)) 
##           0%          10%          20%          30%          40%          50% 
## -42.02999493  -2.87352237   0.00000000   0.04663897   2.46136429   5.19906939 
##          60%          70%          80%          90%         100% 
##   9.34748067  13.67239987  18.11758968  24.02085190 148.15226168

p.breaks <- c(-43,-3,0,3,5,9,13,18,24,148)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")
map_list <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = pr.percent))+
    borders(colour = "grey60")+
   scale_fill_gradient2(high = "#08519c", mid = "white", low = "#b2182b", midpoint = 0, transform = "pseudo_log", breaks = c(-10, -5, -1, 0, 1, 5, 10, 15, 30))+
      labs(fill = "Precipitation anomalies in %", title = i)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = pr.percent))+
      borders(colour = "grey60")+
   scale_fill_gradient2(high = "#08519c", mid = "white", low = "#b2182b", midpoint = 0, transform = "pseudo_log", breaks = c(-10, -5, -1, 0, 1, 5, 10, 15, 30))+
      labs(fill = "Precipitation anomalies in %", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

12.2.2.2 Binned anomalies


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = pr.percent))+
      borders(colour = "grey60")+
   scale_fill_binned_divergingx(mid = 0, palette = "RdBu", rev = F, n_interp = 7, breaks = c(-20, -10, 0, 10, 20, 30, 40, 50, 60))+
      labs(fill = "Precipitation anomalies in %", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

BR <- brewer.pal(n = 11, name ="RdBu")
coltemp <- c("[-20,-10]" = BR[4], "(-10,0]" = BR[6], "(0,10]" = BR[6], "(10,20]" = BR[8], "(20,30]" = BR[9], "(30,40]" = BR[10], "(40,50]" = BR[11], "(50,60]" = "grey20")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot(data = subset(cmip6s, period == j & model == i), 
                aes(x=lon, y = lat,  fill = cut(pr.percent, breaks = c(-20, -10, 0, 10, 20, 30, 40, 50, 60), include.lowest = T))) + 
      geom_tile(show.legend = T)+
      borders(colour = "grey60")+
      scale_fill_manual(values = coltemp, drop = F, na.translate = T)+
      labs(fill = "Precipitation anomalies in %", title = index)+
      theme_void()+ylim(-55,90)+
     theme(legend.position = "bottom", legend.text =  element_text(angle = 0))+
    guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

BR <- brewer.pal(n = 11, name ="RdBu")
coltemp <- c("[-20,-10]" = BR[4], "(-10,0]" = BR[5], "(0,10]" = BR[7], "(10,20]" = BR[8], "(20,30]" = BR[9], "(30,40]" = BR[10], "(40,50]" = BR[11], "(50,60]" = "grey20")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot(data = subset(cmip6s, period == j & model == i), 
                aes(x=lon, y = lat,  fill = cut(pr.percent, breaks = c(-20, -10, 0, 10, 20, 30, 40, 50, 60), include.lowest = T))) + 
      geom_tile(show.legend = T)+
      borders(colour = "grey60")+
      scale_fill_manual(values = coltemp, drop = F, na.translate = T)+
      labs(fill = "Precipitation anomalies in %", title = index)+
      theme_void()+ylim(-55,90)+
     theme(legend.position = "bottom", legend.text =  element_text(angle = 0))+
    guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

13 Seasonnal aridity

In CMIP6, the precipitation (pr) and evapotranspiration (ET0) are in mm/day.

cmip6 <- read.table("cmip6_AI.txt")

cmip6$AI_01 <- with(cmip6, pr_01/ET0_01)
cmip6$AI_02 <- with(cmip6, pr_02/ET0_02)
cmip6$AI_03 <- with(cmip6, pr_03/ET0_03)
cmip6$AI_04 <- with(cmip6, pr_04/ET0_04)
cmip6$AI_05 <- with(cmip6, pr_05/ET0_05)
cmip6$AI_06 <- with(cmip6, pr_06/ET0_06)
cmip6$AI_07 <- with(cmip6, pr_07/ET0_07)
cmip6$AI_08 <- with(cmip6, pr_08/ET0_08)
cmip6$AI_09 <- with(cmip6, pr_09/ET0_09)
cmip6$AI_10 <- with(cmip6, pr_10/ET0_10)
cmip6$AI_11 <- with(cmip6, pr_11/ET0_11)
cmip6$AI_12 <- with(cmip6, pr_12/ET0_12)
cmip6.monthly.ref <- cmip6 %>% filter(period == "1970_2000") %>% 
  select(c("lon","lat","model","period","source", "Continent","Name",
           "AI_01","pr_01","ET0_01",
           "AI_02","pr_02","ET0_02",
           "AI_03","pr_03","ET0_03",
           "AI_04","pr_04","ET0_04",
           "AI_05","pr_05","ET0_05",
           "AI_06","pr_06","ET0_06",
           "AI_07","pr_07","ET0_07",
           "AI_08","pr_08","ET0_08",
           "AI_09","pr_09","ET0_09",
           "AI_10","pr_10","ET0_10",
           "AI_11","pr_11","ET0_11",
           "AI_12","pr_12","ET0_12"
           )) %>% 
  setNames(c("lon","lat","model","period","source","Continent","Name",
          "AI_01.ref","pr_01.ref","ET0_01.ref",
           "AI_02.ref","pr_02.ref","ET0_02.ref",
           "AI_03.ref","pr_03.ref","ET0_03.ref",
           "AI_04.ref","pr_04.ref","ET0_04.ref",
           "AI_05.ref","pr_05.ref","ET0_05.ref",
           "AI_06.ref","pr_06.ref","ET0_06.ref",
           "AI_07.ref","pr_07.ref","ET0_07.ref",
           "AI_08.ref","pr_08.ref","ET0_08.ref",
           "AI_09.ref","pr_09.ref","ET0_09.ref",
           "AI_10.ref","pr_10.ref","ET0_10.ref",
           "AI_11.ref","pr_11.ref","ET0_11.ref",
           "AI_12.ref","pr_12.ref","ET0_12.ref"
           ))

cmip6monthly <- merge(
  select(cmip6, c("lon","lat","source","model","period","source", "Continent","Name",
            "AI_01","pr_01","ET0_01",
           "AI_02","pr_02","ET0_02",
           "AI_03","pr_03","ET0_03",
           "AI_04","pr_04","ET0_04",
           "AI_05","pr_05","ET0_05",
           "AI_06","pr_06","ET0_06",
           "AI_07","pr_07","ET0_07",
           "AI_08","pr_08","ET0_08",
           "AI_09","pr_09","ET0_09",
           "AI_10","pr_10","ET0_10",
           "AI_11","pr_11","ET0_11",
           "AI_12","pr_12","ET0_12")),
  select(cmip6.monthly.ref, c("lon","lat","source",
           "AI_01.ref","pr_01.ref","ET0_01.ref",
           "AI_02.ref","pr_02.ref","ET0_02.ref",
           "AI_03.ref","pr_03.ref","ET0_03.ref",
           "AI_04.ref","pr_04.ref","ET0_04.ref",
           "AI_05.ref","pr_05.ref","ET0_05.ref",
           "AI_06.ref","pr_06.ref","ET0_06.ref",
           "AI_07.ref","pr_07.ref","ET0_07.ref",
           "AI_08.ref","pr_08.ref","ET0_08.ref",
           "AI_09.ref","pr_09.ref","ET0_09.ref",
           "AI_10.ref","pr_10.ref","ET0_10.ref",
           "AI_11.ref","pr_11.ref","ET0_11.ref",
           "AI_12.ref","pr_12.ref","ET0_12.ref")),
  by = c("lon","lat","source"), all.x = T) 

write.table(cmip6monthly,"cmip6monthly.txt")
cmip6monthly <- read.table("cmip6monthly.txt")

cmip6smonthly <- cmip6monthly %>% 
  group_by(lon,lat, Continent, Name, model, period) %>% 
  dplyr::summarise(AI_01_mean = mean(AI_01, na.rm = T), AI_01_sd = sd(AI_01, na.rm = T),
            AI_02_mean = mean(AI_02, na.rm = T), AI_02_sd = sd(AI_02, na.rm = T),
            AI_03_mean = mean(AI_03, na.rm = T), AI_03_sd = sd(AI_03, na.rm = T),
            AI_04_mean = mean(AI_04, na.rm = T), AI_04_sd = sd(AI_04, na.rm = T),
            AI_05_mean = mean(AI_05, na.rm = T), AI_05_sd = sd(AI_05, na.rm = T),
            AI_06_mean = mean(AI_06, na.rm = T), AI_06_sd = sd(AI_06, na.rm = T),
            AI_07_mean = mean(AI_07, na.rm = T), AI_07_sd = sd(AI_07, na.rm = T),
            AI_08_mean = mean(AI_08, na.rm = T), AI_08_sd = sd(AI_08, na.rm = T),
            AI_09_mean = mean(AI_09, na.rm = T), AI_09_sd = sd(AI_09, na.rm = T),
            AI_10_mean = mean(AI_10, na.rm = T), AI_10_sd = sd(AI_10, na.rm = T),
            AI_11_mean = mean(AI_11, na.rm = T), AI_11_sd = sd(AI_11, na.rm = T),
            AI_12_mean = mean(AI_12, na.rm = T), AI_11_sd = sd(AI_12, na.rm = T),
            ET0_01_mean = mean(ET0_01*31, na.rm = T), ET0_01_sd = sd(ET0_01*31, na.rm = T),
            ET0_02_mean = mean(ET0_02*38, na.rm = T), ET0_02_sd = sd(ET0_02*28, na.rm = T),
            ET0_03_mean = mean(ET0_03*31, na.rm = T), ET0_03_sd = sd(ET0_03*31, na.rm = T),
            ET0_04_mean = mean(ET0_04*30, na.rm = T), ET0_04_sd = sd(ET0_04*30, na.rm = T),
            ET0_05_mean = mean(ET0_05*31, na.rm = T), ET0_05_sd = sd(ET0_05*31, na.rm = T),
            ET0_06_mean = mean(ET0_06*30, na.rm = T), ET0_06_sd = sd(ET0_06*30, na.rm = T),
            ET0_07_mean = mean(ET0_07*31, na.rm = T), ET0_07_sd = sd(ET0_07*31, na.rm = T),
            ET0_08_mean = mean(ET0_08*31, na.rm = T), ET0_08_sd = sd(ET0_08*31, na.rm = T),
            ET0_09_mean = mean(ET0_09*30, na.rm = T), ET0_09_sd = sd(ET0_09*30, na.rm = T),
            ET0_10_mean = mean(ET0_10*31, na.rm = T), ET0_10_sd = sd(ET0_10*31, na.rm = T),
            ET0_11_mean = mean(ET0_11*30, na.rm = T), ET0_11_sd = sd(ET0_11*30, na.rm = T),
            ET0_12_mean = mean(ET0_12*31, na.rm = T), ET0_11_sd = sd(ET0_12*31, na.rm = T),
            pr_01_mean = mean(pr_01*31, na.rm = T), pr_01_sd = sd(pr_01*31, na.rm = T),
            pr_02_mean = mean(pr_02*28, na.rm = T), pr_02_sd = sd(pr_02*28, na.rm = T),
            pr_03_mean = mean(pr_03*31, na.rm = T), pr_03_sd = sd(pr_03*31, na.rm = T),
            pr_04_mean = mean(pr_04*30, na.rm = T), pr_04_sd = sd(pr_04*30, na.rm = T),
            pr_05_mean = mean(pr_05*31, na.rm = T), pr_05_sd = sd(pr_05*31, na.rm = T),
            pr_06_mean = mean(pr_06*30, na.rm = T), pr_06_sd = sd(pr_06*30, na.rm = T),
            pr_07_mean = mean(pr_07*31, na.rm = T), pr_07_sd = sd(pr_07*31, na.rm = T),
            pr_08_mean = mean(pr_08*31, na.rm = T), pr_08_sd = sd(pr_08*31, na.rm = T),
            pr_09_mean = mean(pr_09*30, na.rm = T), pr_09_sd = sd(pr_09*30, na.rm = T),
            pr_10_mean = mean(pr_10*31, na.rm = T), pr_10_sd = sd(pr_10*31, na.rm = T),
            pr_11_mean = mean(pr_11*30, na.rm = T), pr_11_sd = sd(pr_11*30, na.rm = T),
            pr_12_mean = mean(pr_12*31, na.rm = T), pr_11_sd = sd(pr_12*31, na.rm = T)
            ) %>%
 ungroup() %>% as.data.frame()

cmip6smonthly.ref <- cmip6smonthly %>% filter(period == "1970_2000") %>% 
  select(c("lon","lat","model","period", "Continent","Name",
           "AI_01_mean","ET0_01_mean","pr_01_mean",
           "AI_02_mean","ET0_02_mean","pr_02_mean",
           "AI_03_mean","ET0_03_mean","pr_03_mean",
           "AI_04_mean","ET0_04_mean","pr_04_mean",
           "AI_05_mean","ET0_05_mean","pr_05_mean",
           "AI_06_mean","ET0_06_mean","pr_06_mean",
           "AI_07_mean","ET0_07_mean","pr_07_mean",
           "AI_08_mean","ET0_08_mean","pr_08_mean",
           "AI_09_mean","ET0_09_mean","pr_09_mean",
           "AI_10_mean","ET0_10_mean","pr_10_mean",
           "AI_11_mean","ET0_11_mean","pr_11_mean",
           "AI_12_mean","ET0_12_mean","pr_12_mean",
           )) %>% 
  setNames(c("lon","lat","model","period","Continent","Name",
          "AI_01_mean.ref","ET0_01_mean.ref","pr_01_mean.ref",
           "AI_02_mean.ref","ET0_02_mean.ref","pr_02_mean.ref",
           "AI_03_mean.ref","ET0_03_mean.ref","pr_03_mean.ref",
           "AI_04_mean.ref","ET0_04_mean.ref","pr_04_mean.ref",
           "AI_05_mean.ref","ET0_05_mean.ref","pr_05_mean.ref",
           "AI_06_mean.ref","ET0_06_mean.ref","pr_06_mean.ref",
           "AI_07_mean.ref","ET0_07_mean.ref","pr_07_mean.ref",
           "AI_08_mean.ref","ET0_08_mean.ref","pr_08_mean.ref",
           "AI_09_mean.ref","ET0_09_mean.ref","pr_09_mean.ref",
           "AI_10_mean.ref","ET0_10_mean.ref","pr_10_mean.ref",
           "AI_11_mean.ref","ET0_11_mean.ref","pr_11_mean.ref",
           "AI_12_mean.ref","ET0_12_mean.ref","pr_12_mean.ref"
           ))

cmip6smonthly <- merge(cmip6smonthly,cmip6smonthly.ref, by = c("lon","lat","model","period", "Continent","Name"), all.x = T)

# CATEGORIES 01
cmip6smonthly$cat.AI_01 <- cmip6smonthly$AI_01_mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_01[which(cmip6smonthly$ET0_01_mean < 400/12)] <- "Cold"

cmip6smonthly$cat.AI_01.ref <- cmip6smonthly$AI_01_mean.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_01.ref[which(cmip6smonthly$ET0_01_mean.ref < 400/12)] <- "Cold"

# CATEGORIES 02
cmip6smonthly$cat.AI_02 <- cmip6smonthly$AI_02_mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_02[which(cmip6smonthly$ET0_02_mean < 400/12)] <- "Cold"

cmip6smonthly$cat.AI_02.ref <- cmip6smonthly$AI_02_mean.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_02.ref[which(cmip6smonthly$ET0_02_mean.ref < 400/12)] <- "Cold"

# CATEGORIES 03
cmip6smonthly$cat.AI_03 <- cmip6smonthly$AI_03_mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_03[which(cmip6smonthly$ET0_03_mean < 400/12)] <- "Cold"

cmip6smonthly$cat.AI_03.ref <- cmip6smonthly$AI_03_mean.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_03.ref[which(cmip6smonthly$ET0_03_mean.ref < 400/12)] <- "Cold"

# CATEGORIES 04
cmip6smonthly$cat.AI_04 <- cmip6smonthly$AI_04_mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_04[which(cmip6smonthly$ET0_04_mean < 400/12)] <- "Cold"

cmip6smonthly$cat.AI_04.ref <- cmip6smonthly$AI_05_mean.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_04.ref[which(cmip6smonthly$ET0_05_mean.ref < 400/12)] <- "Cold"

# CATEGORIES 05
cmip6smonthly$cat.AI_05 <- cmip6smonthly$AI_05_mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_05[which(cmip6smonthly$ET0_05_mean < 400/12)] <- "Cold"

cmip6smonthly$cat.AI_05.ref <- cmip6smonthly$AI_05_mean.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_05.ref[which(cmip6smonthly$ET0_05_mean.ref < 400/12)] <- "Cold"

# CATEGORIES 06
cmip6smonthly$cat.AI_06 <- cmip6smonthly$AI_06_mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_06[which(cmip6smonthly$ET0_06_mean < 400/12)] <- "Cold"

cmip6smonthly$cat.AI_06.ref <- cmip6smonthly$AI_06_mean.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_06.ref[which(cmip6smonthly$ET0_06_mean.ref < 400/12)] <- "Cold"

# CATEGORIES 07
cmip6smonthly$cat.AI_07 <- cmip6smonthly$AI_07_mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_07[which(cmip6smonthly$ET0_07_mean < 400/12)] <- "Cold"

cmip6smonthly$cat.AI_07.ref <- cmip6smonthly$AI_07_mean.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_07.ref[which(cmip6smonthly$ET0_07_mean.ref < 400/12)] <- "Cold"

# CATEGORIES 08
cmip6smonthly$cat.AI_08 <- cmip6smonthly$AI_08_mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_08[which(cmip6smonthly$ET0_08_mean < 400/12)] <- "Cold"

cmip6smonthly$cat.AI_08.ref <- cmip6smonthly$AI_08_mean.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_08.ref[which(cmip6smonthly$ET0_08_mean.ref < 400/12)] <- "Cold"

# CATEGORIES 09
cmip6smonthly$cat.AI_09 <- cmip6smonthly$AI_09_mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_09[which(cmip6smonthly$ET0_09_mean < 400/12)] <- "Cold"

cmip6smonthly$cat.AI_09.ref <- cmip6smonthly$AI_09_mean.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_09.ref[which(cmip6smonthly$ET0_09_mean.ref < 400/12)] <- "Cold"

# CATEGORIES 10
cmip6smonthly$cat.AI_10 <- cmip6smonthly$AI_10_mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_10[which(cmip6smonthly$ET0_10_mean < 400/12)] <- "Cold"

cmip6smonthly$cat.AI_10.ref <- cmip6smonthly$AI_10_mean.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_10.ref[which(cmip6smonthly$ET0_10_mean.ref < 400/12)] <- "Cold"

# CATEGORIES 11
cmip6smonthly$cat.AI_11 <- cmip6smonthly$AI_11_mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_11[which(cmip6smonthly$ET0_11_mean < 400/12)] <- "Cold"

cmip6smonthly$cat.AI_11.ref <- cmip6smonthly$AI_11_mean.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_11.ref[which(cmip6smonthly$ET0_11_mean.ref < 400/12)] <- "Cold"

# CATEGORIES 12
cmip6smonthly$cat.AI_12 <- cmip6smonthly$AI_12_mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_12[which(cmip6smonthly$ET0_12_mean < 400/12)] <- "Cold"

cmip6smonthly$cat.AI_12.ref <- cmip6smonthly$AI_12_mean.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6smonthly$cat.AI_12.ref[which(cmip6smonthly$ET0_12_mean.ref < 400/12)] <- "Cold"

write.table(cmip6smonthly, "cmip6smonthly.txt")
cmip6smonthly <- read.table("cmip6smonthly.txt")

cmip6smonthly$cat.AI_01 <- factor(cmip6smonthly$cat.AI_01, levels = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))

ggarrange(
ggplot() + 
    geom_tile(data = subset(cmip6smonthly, period == "1970_2000" & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI_01))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "January", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom"),

ggplot() + 
    geom_tile(data = subset(cmip6smonthly, period == "1970_2000" & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI_02))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "February", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom"),

ggplot() + 
    geom_tile(data = subset(cmip6smonthly, period == "1970_2000" & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI_03))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "March", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom"),

ggplot() + 
    geom_tile(data = subset(cmip6smonthly, period == "1970_2000" & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI_04))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "April", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom"),

ggplot() + 
    geom_tile(data = subset(cmip6smonthly, period == "1970_2000" & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI_05))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "May", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom"),

ggplot() + 
    geom_tile(data = subset(cmip6smonthly, period == "1970_2000" & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI_06))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "June", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom"),

ggplot() + 
    geom_tile(data = subset(cmip6smonthly, period == "1970_2000" & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI_07))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "July", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom"),

ggplot() + 
    geom_tile(data = subset(cmip6smonthly, period == "1970_2000" & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI_08))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "August", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom"),

ggplot() + 
    geom_tile(data = subset(cmip6smonthly, period == "1970_2000" & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI_09))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "September", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom"),

ggplot() + 
    geom_tile(data = subset(cmip6smonthly, period == "1970_2000" & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI_10))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "October", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom"),

ggplot() + 
    geom_tile(data = subset(cmip6smonthly, period == "1970_2000" & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI_11))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "November", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom"),

ggplot() + 
    geom_tile(data = subset(cmip6smonthly, period == "1970_2000" & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI_12))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = "December", fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom"),

ncol = 3, nrow = 4, common.legend = T, legend = "bottom"
)

13.1 Evolution of monthly aridity by continent

cmip6smm <- cmip6smonthly %>% group_by(model, period, Continent) %>%
  summarise(AI_01 = median(AI_01_mean, na.rm = T), cov_AI_01 = sd(AI_01_mean, na.rm = T)/mean(AI_01_mean, na.rm = T),
            AI_02 = median(AI_02_mean, na.rm = T), cov_AI_02 = sd(AI_02_mean, na.rm = T)/mean(AI_02_mean, na.rm = T),
            AI_03 = median(AI_03_mean, na.rm = T), cov_AI_03 = sd(AI_03_mean, na.rm = T)/mean(AI_03_mean, na.rm = T),
            AI_04 = median(AI_04_mean, na.rm = T), cov_AI_04 = sd(AI_04_mean, na.rm = T)/mean(AI_04_mean, na.rm = T),
            AI_05 = median(AI_05_mean, na.rm = T), cov_AI_05 = sd(AI_05_mean, na.rm = T)/mean(AI_05_mean, na.rm = T),
            AI_06 = median(AI_06_mean, na.rm = T), cov_AI_06 = sd(AI_06_mean, na.rm = T)/mean(AI_06_mean, na.rm = T),
            AI_07 = median(AI_07_mean, na.rm = T), cov_AI_07 = sd(AI_07_mean, na.rm = T)/mean(AI_07_mean, na.rm = T),
            AI_08 = median(AI_08_mean, na.rm = T), cov_AI_08 = sd(AI_08_mean, na.rm = T)/mean(AI_08_mean, na.rm = T),
            AI_09 = median(AI_09_mean, na.rm = T), cov_AI_09 = sd(AI_09_mean, na.rm = T)/mean(AI_09_mean, na.rm = T),
            AI_10 = median(AI_10_mean, na.rm = T), cov_AI_10 = sd(AI_10_mean, na.rm = T)/mean(AI_10_mean, na.rm = T),
            AI_11 = median(AI_11_mean, na.rm = T), cov_AI_11 = sd(AI_11_mean, na.rm = T)/mean(AI_11_mean, na.rm = T),
            AI_12 = median(AI_12_mean, na.rm = T), cov_AI_12 = sd(AI_12_mean, na.rm = T)/mean(AI_12_mean, na.rm = T)) %>%
 ungroup() %>% as.data.frame()

# cmip6smm <- cmip6smonthly %>% group_by(model, period, Continent) %>% 
#   summarise(AI_01 = mean(AI_01_mean, na.rm = T),
#             AI_02 = mean(AI_02_mean, na.rm = T),
#             AI_03 = mean(AI_03_mean, na.rm = T),
#             AI_04 = mean(AI_04_mean, na.rm = T),
#             AI_05 = mean(AI_05_mean, na.rm = T),
#             AI_06 = mean(AI_06_mean, na.rm = T),
#             AI_07 = mean(AI_07_mean, na.rm = T),
#             AI_08 = mean(AI_08_mean, na.rm = T),
#             AI_09 = mean(AI_09_mean, na.rm = T),
#             AI_10 = mean(AI_10_mean, na.rm = T),
#             AI_11 = mean(AI_11_mean, na.rm = T),
#             AI_12 = mean(AI_12_mean, na.rm = T)) %>%
#  ungroup() %>% as.data.frame()

long_AI <- cmip6smm %>% select(c("model","period","Continent","AI_01","AI_02","AI_03","AI_04","AI_05","AI_06","AI_07","AI_08","AI_09","AI_10","AI_11","AI_12")) %>%
                                 tidyr::pivot_longer(
                                   cols = starts_with("AI_"),
                                   names_to= "Month",
                                   values_to = "AI")

long_AI$Month <- revalue(long_AI$Month, c("AI_01" = "January","AI_02" = "February", "AI_03" = "March", "AI_04" = "April","AI_05" ="May","AI_06" = "June", "AI_07" ="July","AI_08" = "August","AI_09" = "September","AI_10" = "October","AI_11" = "November","AI_12" = "December"))
                               
long_sdAI <- cmip6smm %>% 
  select(c("model","period","Continent","cov_AI_01","cov_AI_02","cov_AI_03","cov_AI_04","cov_AI_05","cov_AI_06","cov_AI_07","cov_AI_08","cov_AI_09","cov_AI_10","cov_AI_11","cov_AI_12")) %>%
  tidyr::pivot_longer(
  cols = starts_with("cov_AI_"),
  names_to= "Month",
  values_to = "cov_AI"
)

long_sdAI$Month <- revalue(long_sdAI$Month, c("cov_AI_01" = "January","cov_AI_02" = "February", "cov_AI_03" = "March", "cov_AI_04" = "April","cov_AI_05" ="May","cov_AI_06" = "June", "cov_AI_07" ="July","cov_AI_08" = "August","cov_AI_09" = "September","cov_AI_10" = "October","cov_AI_11" = "November","cov_AI_12" = "December"))

cmip6sml <- merge(long_AI,long_sdAI, by = c("model","period","Continent","Month"))

cmip6sml$Month <- factor(cmip6sml$Month, c("January","February","March","April","May","June","July","August","September","October","November","December"))

col.model <- c("historical" = "#4575B4", "SSP245" = "#FEE08B", "SSP370" = "#FDAE61", "SSP585" = "#D73027")

cmip6sml$NMonth <- cmip6sml$Month %>% revalue(c("January" = 1, "February" =2, "March"= 3, "April"= 4,"May" = 5, "June" = 6,"July" = 7, "August"= 8,"September" = 9, "October" = 10,"November" = 11, "December" = 12))

for(i in  c("NORTH-AMERICA","CENTRAL-AMERICA","SOUTH-AMERICA","EUROPE","EUROPE-AFRICA","AFRICA","ASIA","OCEANIA")){
g <- ggplot(subset(cmip6sml, Continent == i & period %in% c("1970_2000","2070_2100")))+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0, ymax = 0.03), fill = "#D73027", alpha = 0.02)+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0.03, ymax = 0.2), fill = "#FDAE61", alpha = 0.02)+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0.2, ymax = 0.5), fill = "#FEE08B", alpha = 0.02)+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0.5, ymax = 0.65), fill = "#ABD9E9", alpha = 0.02)+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0.65, ymax = Inf), fill = "#4575B4", alpha = 0.02)+
  geom_point(aes(x=Month, y = AI, size = cov_AI, col = model), alpha = 0.5)+
  geom_line(aes(x=Month, y = AI, col = model, group = model, linetype = model), alpha = 0.5)+
  scale_color_manual(values = col.model)+
  scale_size_continuous(transform = "pseudo_log")+
  labs(y = "Median aridity index", size = "Coefficient of variance", col = "Model", linetype = "Model", title = i)+
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
print(g)
}

13.2 Evolution of monthly evapotranspiration by continent

cmip6smm_et0 <- cmip6smonthly %>% group_by(model, period, Continent) %>%
  summarise(ET0_01 = median(ET0_01_mean, na.rm = T), cov_ET0_01 = sd(ET0_01_mean, na.rm = T)/mean(ET0_01_mean, na.rm = T),
            ET0_02 = median(ET0_02_mean, na.rm = T), cov_ET0_02 = sd(ET0_02_mean, na.rm = T)/mean(ET0_02_mean, na.rm = T),
            ET0_03 = median(ET0_03_mean, na.rm = T), cov_ET0_03 = sd(ET0_03_mean, na.rm = T)/mean(ET0_03_mean, na.rm = T),
            ET0_04 = median(ET0_04_mean, na.rm = T), cov_ET0_04 = sd(ET0_04_mean, na.rm = T)/mean(ET0_04_mean, na.rm = T),
            ET0_05 = median(ET0_05_mean, na.rm = T), cov_ET0_05 = sd(ET0_05_mean, na.rm = T)/mean(ET0_05_mean, na.rm = T),
            ET0_06 = median(ET0_06_mean, na.rm = T), cov_ET0_06 = sd(ET0_06_mean, na.rm = T)/mean(ET0_06_mean, na.rm = T),
            ET0_07 = median(ET0_07_mean, na.rm = T), cov_ET0_07 = sd(ET0_07_mean, na.rm = T)/mean(ET0_07_mean, na.rm = T),
            ET0_08 = median(ET0_08_mean, na.rm = T), cov_ET0_08 = sd(ET0_08_mean, na.rm = T)/mean(ET0_08_mean, na.rm = T),
            ET0_09 = median(ET0_09_mean, na.rm = T), cov_ET0_09 = sd(ET0_09_mean, na.rm = T)/mean(ET0_09_mean, na.rm = T),
            ET0_10 = median(ET0_10_mean, na.rm = T), cov_ET0_10 = sd(ET0_10_mean, na.rm = T)/mean(ET0_10_mean, na.rm = T),
            ET0_11 = median(ET0_11_mean, na.rm = T), cov_ET0_11 = sd(ET0_11_mean, na.rm = T)/mean(ET0_11_mean, na.rm = T),
            ET0_12 = median(ET0_12_mean, na.rm = T), cov_ET0_12 = sd(ET0_12_mean, na.rm = T)/mean(ET0_12_mean, na.rm = T)) %>%
 ungroup() %>% as.data.frame()

# cmip6smm <- cmip6smonthly %>% group_by(model, period, Continent) %>% 
#   summarise(ET0_01 = mean(ET0_01_mean, na.rm = T),
#             ET0_02 = mean(ET0_02_mean, na.rm = T),
#             ET0_03 = mean(ET0_03_mean, na.rm = T),
#             ET0_04 = mean(ET0_04_mean, na.rm = T),
#             ET0_05 = mean(ET0_05_mean, na.rm = T),
#             ET0_06 = mean(ET0_06_mean, na.rm = T),
#             ET0_07 = mean(ET0_07_mean, na.rm = T),
#             ET0_08 = mean(ET0_08_mean, na.rm = T),
#             ET0_09 = mean(ET0_09_mean, na.rm = T),
#             ET0_10 = mean(ET0_10_mean, na.rm = T),
#             ET0_11 = mean(ET0_11_mean, na.rm = T),
#             ET0_12 = mean(ET0_12_mean, na.rm = T)) %>%
#  ungroup() %>% as.data.frame()

long_ET0 <- cmip6smm_et0 %>% select(c("model","period","Continent","ET0_01","ET0_02","ET0_03","ET0_04","ET0_05","ET0_06","ET0_07","ET0_08","ET0_09","ET0_10","ET0_11","ET0_12")) %>%
                                 tidyr::pivot_longer(
                                   cols = starts_with("ET0_"),
                                   names_to= "Month",
                                   values_to = "ET0")

long_ET0$Month <- revalue(long_AI$Month, c("ET0_01" = "January","ET0_02" = "February", "ET0_03" = "March", "ET0_04" = "April","ET0_05" ="May","ET0_06" = "June", "ET0_07" ="July","ET0_08" = "August","ET0_09" = "September","ET0_10" = "October","ET0_11" = "November","ET0_12" = "December"))
                               
long_sdET0 <- cmip6smm_et0 %>% 
  select(c("model","period","Continent","cov_ET0_01","cov_ET0_02","cov_ET0_03","cov_ET0_04","cov_ET0_05","cov_ET0_06","cov_ET0_07","cov_ET0_08","cov_ET0_09","cov_ET0_10","cov_ET0_11","cov_ET0_12")) %>%
  tidyr::pivot_longer(
  cols = starts_with("cov_ET0_"),
  names_to= "Month",
  values_to = "cov_ET0"
)

long_sdET0$Month <- revalue(long_sdET0$Month, c("cov_ET0_01" = "January","cov_ET0_02" = "February", "cov_ET0_03" = "March", "cov_ET0_04" = "April","cov_ET0_05" ="May","cov_ET0_06" = "June", "cov_ET0_07" ="July","cov_ET0_08" = "August","cov_ET0_09" = "September","cov_ET0_10" = "October","cov_ET0_11" = "November","cov_ET0_12" = "December"))

cmip6sml_et0 <- merge(long_ET0,long_sdET0, by = c("model","period","Continent","Month"))

cmip6sml_et0$Month <- factor(cmip6sml_et0$Month, c("January","February","March","April","May","June","July","August","September","October","November","December"))

col.model <- c("historical" = "#4575B4", "SSP245" = "#FEE08B", "SSP370" = "#FDAE61", "SSP585" = "#D73027")

cmip6sml_et0$NMonth <- cmip6sml_et0$Month %>% revalue(c("January" = 1, "February" =2, "March"= 3, "April"= 4,"May" = 5, "June" = 6,"July" = 7, "August"= 8,"September" = 9, "October" = 10,"November" = 11, "December" = 12))

for(i in  c("NORTH-AMERICA","CENTRAL-AMERICA","SOUTH-AMERICA","EUROPE","EUROPE-AFRICA","AFRICA","ASIA","OCEANIA")){
g <- ggplot(subset(cmip6sml_et0, Continent == i & period %in% c("1970_2000","2070_2100")))+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0, ymax = 0.03), fill = "#D73027", alpha = 0.02)+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0.03, ymax = 0.2), fill = "#FDAE61", alpha = 0.02)+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0.2, ymax = 0.5), fill = "#FEE08B", alpha = 0.02)+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0.5, ymax = 0.65), fill = "#ABD9E9", alpha = 0.02)+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0.65, ymax = Inf), fill = "#4575B4", alpha = 0.02)+
  geom_point(aes(x=Month, y = ET0, size = cov_ET0, col = model), alpha = 0.5)+
  geom_line(aes(x=Month, y = ET0, col = model, group = model, linetype = model), alpha = 0.5)+
  scale_color_manual(values = col.model)+
  scale_size_continuous(transform = "pseudo_log")+
  labs(y = "Median monthly evapotranspiration, mm", size = "Coefficient of variance", col = "Model", linetype = "Model", title = i)+
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
print(g)
}

13.3 Evolution of monthly precipitation by continent

cmip6smm_pr <- cmip6smonthly %>% group_by(model, period, Continent) %>%
  summarise(pr_01 = median(pr_01_mean, na.rm = T), cov_pr_01 = sd(pr_01_mean, na.rm = T)/mean(pr_01_mean, na.rm = T),
            pr_02 = median(pr_02_mean, na.rm = T), cov_pr_02 = sd(pr_02_mean, na.rm = T)/mean(pr_02_mean, na.rm = T),
            pr_03 = median(pr_03_mean, na.rm = T), cov_pr_03 = sd(pr_03_mean, na.rm = T)/mean(pr_03_mean, na.rm = T),
            pr_04 = median(pr_04_mean, na.rm = T), cov_pr_04 = sd(pr_04_mean, na.rm = T)/mean(pr_04_mean, na.rm = T),
            pr_05 = median(pr_05_mean, na.rm = T), cov_pr_05 = sd(pr_05_mean, na.rm = T)/mean(pr_05_mean, na.rm = T),
            pr_06 = median(pr_06_mean, na.rm = T), cov_pr_06 = sd(pr_06_mean, na.rm = T)/mean(pr_06_mean, na.rm = T),
            pr_07 = median(pr_07_mean, na.rm = T), cov_pr_07 = sd(pr_07_mean, na.rm = T)/mean(pr_07_mean, na.rm = T),
            pr_08 = median(pr_08_mean, na.rm = T), cov_pr_08 = sd(pr_08_mean, na.rm = T)/mean(pr_08_mean, na.rm = T),
            pr_09 = median(pr_09_mean, na.rm = T), cov_pr_09 = sd(pr_09_mean, na.rm = T)/mean(pr_09_mean, na.rm = T),
            pr_10 = median(pr_10_mean, na.rm = T), cov_pr_10 = sd(pr_10_mean, na.rm = T)/mean(pr_10_mean, na.rm = T),
            pr_11 = median(pr_11_mean, na.rm = T), cov_pr_11 = sd(pr_11_mean, na.rm = T)/mean(pr_11_mean, na.rm = T),
            pr_12 = median(pr_12_mean, na.rm = T), cov_pr_12 = sd(pr_12_mean, na.rm = T)/mean(pr_12_mean, na.rm = T)) %>%
 ungroup() %>% as.data.frame()

# cmip6smm <- cmip6smonthly %>% group_by(model, period, Continent) %>% 
#   summarise(pr_01 = mean(pr_01_mean, na.rm = T),
#             pr_02 = mean(pr_02_mean, na.rm = T),
#             pr_03 = mean(pr_03_mean, na.rm = T),
#             pr_04 = mean(pr_04_mean, na.rm = T),
#             pr_05 = mean(pr_05_mean, na.rm = T),
#             pr_06 = mean(pr_06_mean, na.rm = T),
#             pr_07 = mean(pr_07_mean, na.rm = T),
#             pr_08 = mean(pr_08_mean, na.rm = T),
#             pr_09 = mean(pr_09_mean, na.rm = T),
#             pr_10 = mean(pr_10_mean, na.rm = T),
#             pr_11 = mean(pr_11_mean, na.rm = T),
#             pr_12 = mean(pr_12_mean, na.rm = T)) %>%
#  ungroup() %>% as.data.frame()

long_pr <- cmip6smm_pr %>% select(c("model","period","Continent","pr_01","pr_02","pr_03","pr_04","pr_05","pr_06","pr_07","pr_08","pr_09","pr_10","pr_11","pr_12")) %>%
                                 tidyr::pivot_longer(
                                   cols = starts_with("pr_"),
                                   names_to= "Month",
                                   values_to = "pr")

long_pr$Month <- revalue(long_AI$Month, c("pr_01" = "January","pr_02" = "February", "pr_03" = "March", "pr_04" = "Apr_il","pr_05" ="May","pr_06" = "June", "pr_07" ="July","pr_08" = "August","pr_09" = "September","pr_10" = "October","pr_11" = "November","pr_12" = "December"))
                               
long_sdpr <- cmip6smm_pr %>% 
  select(c("model","period","Continent","cov_pr_01","cov_pr_02","cov_pr_03","cov_pr_04","cov_pr_05","cov_pr_06","cov_pr_07","cov_pr_08","cov_pr_09","cov_pr_10","cov_pr_11","cov_pr_12")) %>%
  tidyr::pivot_longer(
  cols = starts_with("cov_pr_"),
  names_to= "Month",
  values_to = "cov_pr"
)

long_sdpr$Month <- revalue(long_pr$Month, c("cov_pr_01" = "January","cov_pr_02" = "February", "cov_pr_03" = "March", "cov_pr_04" = "April","cov_pr_05" ="May","cov_pr_06" = "June", "cov_pr_07" ="July","cov_pr_08" = "August","cov_pr_09" = "September","cov_pr_10" = "October","cov_pr_11" = "November","cov_pr_12" = "December"))

cmip6sml_pr <- merge(long_pr,long_sdpr, by = c("model","period","Continent","Month"))

cmip6sml_pr$Month <- factor(cmip6sml_pr$Month, c("January","February","March","April","May","June","July","August","September","October","November","December"))

col.model <- c("historical" = "#4575B4", "SSP245" = "#FEE08B", "SSP370" = "#FDAE61", "SSP585" = "#D73027")

cmip6sml_pr$NMonth <- cmip6sml_pr$Month %>% revalue(c("January" = 1, "February" =2, "March"= 3, "April"= 4,"May" = 5, "June" = 6,"July" = 7, "August"= 8,"September" = 9, "October" = 10,"November" = 11, "December" = 12))

for(i in  c("NORTH-AMERICA","CENTRAL-AMERICA","SOUTH-AMERICA","EUROPE","EUROPE-AFRICA","AFRICA","ASIA","OCEANIA")){
g <- ggplot(subset(cmip6sml_pr, Continent == i & period %in% c("1970_2000","2070_2100")))+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0, ymax = 0.03), fill = "#D73027", alpha = 0.02)+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0.03, ymax = 0.2), fill = "#FDAE61", alpha = 0.02)+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0.2, ymax = 0.5), fill = "#FEE08B", alpha = 0.02)+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0.5, ymax = 0.65), fill = "#ABD9E9", alpha = 0.02)+
  # geom_rect(aes(xmin = 1, xmax = 12, ymin= 0.65, ymax = Inf), fill = "#4575B4", alpha = 0.02)+
  geom_point(aes(x=Month, y = pr, size = cov_pr, col = model), alpha = 0.5)+
  geom_line(aes(x=Month, y = pr, col = model, group = model, linetype = model), alpha = 0.5)+
  scale_color_manual(values = col.model)+
  scale_size_continuous(transform = "pseudo_log")+
  labs(y = "Median monthly precipitation, mm", size = "Coefficient of variance", col = "Model", linetype = "Model", title = i)+
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
print(g)
}